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The QCDNUM program numerically solves the evolution equations for parton den- 

-^ i sities and fragmentation functions in perturbative QCD. Un-polarised parton den- 

^ ■ sities can be evolved up to next-to-next-to-leading order in powers of the strong 

coupling constant, while polarised densities or fragmentation functions can be 
evolved up to next-to-leading order. Other types of evolution can be accessed by 
feeding alternative sets of evolution kernels into the program. A versatile con- 

ly^ ■ volution engine provides tools to compute parton luminosities, cross-sections in 

^^ i hadron-hadron scattering, and deep inelastic structure functions in the zero-mass 

scheme or in generalised mass schemes. Input to these calculations are either the 
QCDNUM evolved densities, or those read in from an external parton density repos- 
itory. Included in the software distribution are packages to calculate zero-mass 

rS ' structure functions in un-polarised deep inelastic scattering, and heavy flavour 

C^ ' contributions to these structure functions in the fixed flavour number scheme. 
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1 Introduction 



In perturbative quantum chromo dynamics (pQCD), a hard hadron-hadron scattering 
cross section is calculated as the convolution of a partonic cross section with the mo- 
mentum distributions of the partons inside the colliding hadrons. These parton distri- 
butions depend on the Bjorken-x variable (fractional momentum of the partons inside 
the hadron) and on a scale /i^ characteristic of the hard scattering process. Whereas 
the x-dependence of the parton densities is non-perturbative, the //^ dependence can be 
described in pQCD by the DGLAP evolution equations pQ. The perturbative expansion 
of the splitting functions in these equations has recently been calculated up to next-to- 
next-to- leading order (NNLO) in powers of the strong coupling constant Os [21 13]. 

QcDNUM is a FORTRAN program that numerically solves the DGLAP evolution equa- 
tions on a discrete grid in x and iJ? . Input to the evolution are the x-dependence of the 
parton densities at some input mass factorisation scale, and an input value of as at some 
input renormalisation scale. To study the scale uncertainties, the renormalisation scale 
can be varied with respect to the mass factorisation scale. All calculations in qcdnum 
are performed in the MS scheme. 

The program was originally developed in 1988 by members of the BCDMS collabora- 
tion [4j for a next-to-leading order (NLO) pQCD analysis of the SLAG and BGDMS 
structure function data ^ . This code was adapted by the NMG for use at low x [6] . A 
complete revision led to the version 16.12 which was used in the QGD fits by ZEUS [7], 
and in a global QGD analysis of deep inelastic scattering data by the present author [8]. 

QcDNUMlT is the NNLO upgrade of qcdnum16. A new evolution algorithm, based 
on quadratic spline interpolation, yields large gains in accuracy and speed; on a 2 GHz 
processor it takes less than 10 ms to evolve over a large kinematic range the full set of 
parton densities at NNLO in the variable flavour number scheme. Qcdnum17 can evolve 
un-polarised parton densities up to NNLO, and polarised densities or fragmentation 
functions up to NLO. Alternative sets of evolution kernels can be fed into the program 
to perform other types of evolution. It is also possible to read a parton density set from 
an external library, instead of evolving these from the input scale. 

A versatile set of convolution routines is provided that can be used to calculate hadron- 
hadron scattering cross-sections, parton luminosities, or deep inelastic structure func- 
tions in either the zero-mass or in generalised mass schemes. Included in the software 
distribution are the zmstf and hqstf add-on packages to compute un-polarised zero- 
mass structure functions and, in the fixed flavour number scheme, the contribution from 
heavy quarks to these structure functions. 

This write-up is organised as follows. In Section [2] we summarise the formalism un- 
derlying the DGLAP evolution of parton densities. The qcdnum numerical method 
is described in Section [31 Details about the program itself and the description of an 
example job can be found in Section HI A subroutine-by-subroutine manual is given in 
Section O The qcdnum convolution engine is presented in Section O The zmstf and 
HQSTF packages are described in the Appendices O and [Dl respectively. 



2 QCD Evolution 

In pQCD, the strong coupling constant Og evolves on the renormalisation scale fi^. The 
starting value is specified at some input scale, which usually is taken to be m|. 

The parton density functions (pdf) evolve on the factorisation scale /i|. The starting 
point of a pdf evolution is given by the x dependence of the pdf at some initial scale /ip. 
The coupled evolution equations that are obeyed by the gluon and the quark densities 
can, to a large extent, be decoupled by writing them in terms of the singlet quark density 
(sum of all active quarks and anti-quarks) and non-singlet densities (orthogonal to the 
singlet in fiavour space). A nice feature of qcdnum is that it automatically takes care 
of the singlet /non-singlet decomposition of a set of pdfs. 

Another input to the QCD evolution is the number of active fiavours Uf which specifies 
how many quark species (d, u, s, . . .) are participating in the QCD dynamics. In the 
fixed flavour number scheme (ffns), Uf is kept fixed throughout the evolution. Input to 
an FFNS evolution are then the gluon density and 2nf (anti-) quark densities at the input 
scale y^Q. In the variable flavour number scheme (VFNS), the flavour thresholds /i^bt ^^^ 
introduced and 3 light quark densities (d, u, s) are, together with the corresponding 
anti-quark densities, specified below the charm threshold /i^. The heavy quarks and 
anti-quarks (c, b, t) are dynamicafly generated by the QCD evolution equations at and 
above their thresholds. Both the ffns and the VFNS are supported by qcdnum. 

The QCD evolution formalism is relatively simple when the renormalisation and fac- 
torisation scales are equal, but it becomes more complicated when /i^ 7^ /i|. Qcdnum 
supports a linear relationship between the two scales. 

In the following sections we describe the evolution of Og and the pdfs, the renormalisation 
scale dependence, the singlet /non-singlet decomposition, and the fiavour schemes. 



2.1 Evolution of the Strong Coupling Constant 

The evolution of the strong coupling constant reads, up to NNLO, 

^ = -X:fta«(A=). (2.1) 

Here /i^ = /i^ is the renormalisation scale and Og = a^j^T^. The /9-functions in (12. ip 
depend on the number Uf of active quarks with pole mass m < fi. In the MS scheme 
they are given by [9l [10] 

51 19 

Pi = nf 

2857 5033 325 2 

02 = n/- H ni. (2.2) 

^^ 16 144 -^ 432 ^ ^ ' 

The leading order (LO) analytical solution of (12. ip can be written as 

^ ^ + /3o In f41 = /3o In f^ V (2.3) 



as(/i2) a^inl) \nlJ VA2 
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In fl2.3p . the parameter A is defined as the scale where the first term on the right-hand 
side vanishes, that is, the scale where as becomes infinite. Beyond LO, the definition of a 
scale parameter is ambiguous so that it is more convenient to take as(m|) as a reference. 
The value of as at any other scale is then obtained from a numerical integration of fl2.ip j^l 
instead of from approximate analytical solutions parametrised in terms of A. 

In the evolution of as, the number of active flavours is set to Uf = 3 below the charm 
threshold /i| = fi^ and is changed from nj to nj + l at the flavour thresholds /i^ = /x^ b t- 
At NNLO, and sometimes also at NLO, there are small discontinuities in the as evolution 



at the flavour thresholds [TT] ; see Section 12.51 for details. 



In Figure [H we plot the evolution of as calculated at LO, NLO and NNLOo Because 
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Figure 1: The strong coupling constant as{fi^) evolved downward from as(m|) = 0.118 in LO (dotted 
curve), NLO (dashed curve) and NNLO (full curve). The inset shows an enlarged view of the NNLO 
discontinuity in as at the charm threshold /i^. 

pQCD breaks down when as becomes large, qcdnum will issue a fatal error when as(/i^) 
exceeds a pre-set limit. For a given value of as(m|), it is clear from the figure that such 
a limit will correspond to larger values of fi"^ at larger perturbative order. 



2.2 The DGLAP Evolution Equations 

The DGLAP evolution equations can be written as 



dlnfi'^ 



E 



J=<1,<1,9 



''^"p^A-^^^')fA^^^^') 



z \z 



(2.4) 



^I thank A. Vogt for providing his 4**^ order Runge-Kutta routine to integrate (|2.ip up to NNLO. 
2With the settings as(m|) = 0.118 and ^c,b,t = (1-5,5, 188) GcV. 



where fi denotes an un-polarised parton number density, Pij are the QCD sphtting 
functions, x is the Bjorken scahng variable and /i^ = /i| is the mass factorisation scale, 
which we assume here to be equal to the renormalisation scale /i|^. The indices i and j 
in ( 12 ■4p run over the parton species i.e., the gluon and Uf active flavours of quarks and 
anti-quarks. In the quark parton model, and also in LO pQCD, the parton densities 
are defined such that f{x,^^)dx is, at a given /z^, the number of partons which carry a 
fraction of the nucleon momentum between x and x + da;. The distribution xf{x, yU^) is 
then the parton momentum densityo Beyond LO there is no such intuitive interpreta- 
tion. The definition of / then depends on the renormalisation and factorisation scheme 
in which the calculations are carried out (MS in qcdnum)o 

Introducing a short-hand notation for the Mellin convolution, 

,;«,](.) -/'^/(f).W-£^/W«0. (2.5) 

we can write (12. 4p in compact form as (we drop the arguments x and //^ in the following) 

j=<i,q,a 

If the X dependencies of the parton densities are known at some scale /ig, they can be 
evolved to other values of /x^ by solving this set of 2nf + 1 coupled integro-differential 
equations. Fortunately, (12. 6p can be considerably simplified by taking the symmetries 
in the splitting functions into account |9]: 

P = P - = P 

P = P- = P 

-' q«g -* qig 2r) ^^ 

p p X pv I ps 

-"^q^qfe ~ -"^q^qfe ~ '^ik-'^qq ~r -r qq 

-* qiqfc ~ -* qiqfc ~ ^ikPqq + Pqq- ['^■'J 

Inserting (12. 7p in (12.60 . we find after some algebra that the singlet quark density 

Uf 

gs = 5^(g. + g^) (2.8) 

obeys an evolution equation coupled to the gluon density 

9 /9s^ _ fPqq Pqi,\ ^ f qs\ ^ ^2.9) 



with Pqq given by 



91n/i2 \gj \^Pgq PgJ \^^ 



^qq = ^<; + Pqq + %(^qq + P^ ' (2-10) 



■^In this section we use the number densities f{x,iJ?). In QCDNUM itself, however, we use xf{x,iJ,'^). 
^In the DIS scheme / is defined such that the LO (quark-parton model) expression for the F2 
structure function is preserved at NLO. But this is true only for F2 and not for Fl and xF^. 
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Likewise, we find that the non-singlet combinations 



?ij = (* ± ?0 - (9j ± 9i) and qv = ^{qi-qi) (2.11' 



i=l 



evolve independently from the gluon and from each other according to 



., 2 P±®^r, and _^ = P,®g,, (2.12) 

(7 In /i^ ■' am /x^ 

with splitting functions defined by 

P^ = P^^^±P^^ and P^ = P^^-P^^^ + nf{P^^-P^^). (2.13) 

The evolution of the g^ is linear in the densities, so that any linear combination of the 
q^j or q^j also evolves according to (I2.12p . 

The splitting functions can be expanded in a perturbative series in Og which presently 
is known up to NNLO. For the four splitting functions Pij in 02. 9p we may write 

P,,(x,/i^) = a3(/.2) P^\x) + ajif,') P^\x) + alif,') P^{x) + 0(af) (2.14) 

where we have set, as in the previous section, Og = as/27r. Note the separation in the 
variables x and /i^ on the right-hand side of (I2.14p . We drop again the arguments x 
and /i^ and write the expansion of the non-singlet splitting functions as 



Pv = a,Pl^^+alP^}^ +alPl^'^ + 0{at). (2.15) 



Pi = a,P^^+alP^2^+alPi^ + 0{at) 

• —1— ( / I 

qq 



Truncating the right-hand side to the appropriate order in Og, it is seen that at LO 
the three types of non-singlet obey the same evolution equations. At NLO, q~j and q^ 
evolve in the same way but different from q^-. At NNLO, all three non-singlets evolve 
differently. 

It is evident from (12.71) . (12.101) and (I2.13P that several splitting functions depend on 
the number of active flavours rtf. This number is set to 3 below /i| = ^1 and changed 
to Uf = (4,5,6) at and above the thresholds /i| = /icbt- ^^ ^ase /i| 7^ /i|^, qcdnum 
adjusts the /i^ thresholds such that Uf changes in both the splitting and the beta 
functions when crossing a threshold; see also Section 12.51 

The LO splitting functions are given in Appendix |Xl Those at NLO can be found in [12] 
(non-singlet) and [13] (singlet) The NNLO splitting functions and their parametrisa- 
tions are given in [2] (non-singlet) and [3] (singlet). The DGLAP equations also apply 
to polarised parton densities and to fragmentation functions (time-like evolution), each 
with their own set of evolution kernels. For the polarised splitting functions up to NLO 
we refer to [2] , and references therein. The time-like evolution of fragmentation func- 
tions at LO is described in [15] , see also Appendix \M The NLO time-like splitting 
functions can be found in |12| and [13]. 



^Two well-known misprints in |13| are: (i) the lower integration limit in the definition of 6*2(0;) must 
read a;/(l+x); (ii) in the expression for ipp the term (10 — 18a; — -ya;^) must read (—10— 18a; — ^x^). 



2.3 Renormalisation Scale Dependence 

In the previous section, we have assumed that the factorisation and renormahsation 
scales are equal. For //| 7^ fi^ we expand as in a Taylor series on a logarithmic scale 
around /x^ 

asifil) = a^if^l) + a',{fil)Ln + - a'^ifi^Ll +... (2.16) 

with Lr = ln(/i|//i|). Using (12. ip to calculate the derivatives in (12.161) . we obtain 

ali^il) = al{^,l) + 0{at). (2.17) 

To calculate the renormalisation scale dependence of the evolved parton densities, the 
powers of Og in the splitting function expansions (I2.14p and (I2.15P are replaced by the 
expressions on the right-hand side of (I2.17p . with the understanding that these are 
truncated to order Og when we evolve at LO, to order Og when we evolve at NLO, and 
to order a^ when we evolve at NNLO. 



2.4 Decomposition into Singlet and Non-singlets 

In this section we describe the transformations between a flavour basis and a singlet/non- 
singlet basis, as is implemented in qcdnum. For this purpose we write an arbitrary 
linear combination of quark and anti-quark densities as 

Uf 

\p) = J2ia,\q,)+(3M), (2.18) 

where the index i runs over the number of active flavours. To make a clear distinction 
between a coefficient and a pdf, we introduce here the ket notation |/) for /(a;, /i^). 

Because a linear combination of non-singlets is again a non-singlet, it follows directly 
from the definition (12.111) that the coefficients of any non-singlet satisfy the constraint 

"/ 
5^(«, + A) = 0, (2.19) 

that is, a non-singlet is — by definition — orthogonal to the singlet in fiavour space. 
It is convenient to define Iq^ ) = \qi) ± \qi) and write the linear combination (12.180 as 

\p) = Y.^K\Qt) + KW))- (2.20) 

The coefficients bf, at and /3j are related by 

bt = ^^, a^ = K + K. P. = K-K- (2.21) 
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We define a basis of singlet, valence, and 2{nf — 1) additional non-singlets by 



i-l 



|e^) = |gs), |ei) = |g^). 



\^f) = E l^f ) - (^ - 1) 1^^^) fo^ 2 < z < rif. (2.22) 



In matrix notation, this transformation can be written as 

\e^) = U\q^), 
where U is the nf x nf sub-matrix of the 6x6 transformation matrix 

/I 1 1 1 1 1 \ 
1-10 
11-2000 
111-300 
1 1 1 1-4 

\ 1 1 1 1 1 -5 / 



(2.23) 



U 



(2.24) 



It is seen that the second to sixth row of (I2.24p are orthogonal to the first row (singlet), 
so that they indeed represent non-singlets as defined by (12.191) . In fact, all rows of U are 
orthogonal to each other, so that scaling by the row-wise norm yields a rotation matrix, 
which has the transpose as its inverse. By scaling back this inverse we obtain 



U-' = U^S' 



Tc2 



(2.25) 



where U is the transpose of U and S is the square of the diagonal scaling matrix: 

-1 

fc=l / 



V 



^ij ~ "ij 



Y."-- 



for i = 1 

5ij/i{i — 1) for i > 1. 



(2.26) 



Using (I2.25P and (I2.26P to invert any rif x nf sub-matrix of (12.240 . it is straight forward 
to show by explicit calculation that 



Ul 



l/uf for j = 1 

-l/j foij = iy^l 

l/j(j-l) forj>z 

otherwise. 



(2.27) 



The inverse of the transformation (12.221) is thus given by 

+E 



lit) 



\ef) , ^ \ef) 



^f 



i=2 
i 



jU - 1) 



+ E V 



ef) 



i> 1. 



We can now write the linear combination \p) on the |e^) basis as 
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(2.28) 



(2.29) 



i=l 



where the coefficients df are related to the bf of (I2.2()p by 

Uf Uf 

df = E ^UJ^^ ^f = E 4^^^- (2-30) 

i=i i=i 

Let the starting values of the DGLAP evolutions be given by the gluon density and 
2n/ arbitrary quark densities, that is, by 2n/ + 1 functions of x at some input scale /Iq. 
We can arrange the input quark densities in a 2n/-diniensional vector |p). Likewise, we 
store the densities \q^) in a vector |q), the \ef) in a vector |e) and the 6^ coefficients 
of each input density in the rows of a 2nj x 2nf matrix B. The flavour decomposition 
of the input densities can then be written as \p) = B\q) and the singlet/non-singlet 
decomposition as 

\p) = BT-'\e) with r=(o ^Y (2.31) 

Provided that B~^ exists (i.e. the input densities are hnearly independent), the starting 
values of the singlet and non-singlet densities are calculated from the inverse relation 

\e)=TB-'^\p). (2.32) 



2.5 Flavour Number Schemes 

QcDNUM supports two evolution schemes, known as the fixed fiavour number scheme 
(ffns) and the variable flavour number scheme (VFNs). 

In the FFNS we assume that Uf quark flavours have zero mass, while those of the remain- 
ing flavours are taken to be inflnitely large. In this way, only Uf flavours participate 
in the QCD dynamics so that in the ffns the value of Uf is simply kept constant for 
all /i^, with 3 < Uf < 6. In the ffns, the input scale /Xq can be chosen anywhere within 
the boundaries of the evolution grid, although one should be careful with backward 
evolution in qcdnum; see Section [X^ 

In the VFNS, the number of flavours changes from nj to Uf + 1 when the factorisation 
scale is equal to the pole mass of the heavy quarks /i^ = m\, h = (c, b,t). A heavy 
quark h is thus considered to be inflnitely massive below fi^ and mass- less above /x^. As 
a consequence, the heavy flavour distributions are zero below their respective thresholds 
and are dynamically generated by the QCD evolution equations at and above /z^. Such 
an abrupt turn-on at a flxed scale is of course unphysical but this poses no problem 
since the parton densities themselves are not observables. The VFNS or ffns parton 
densities evolved with qcdnum are, in fact, valid input to structure function and cross 
section calculations that include mass terms and obey the kinematics of heavy quark 
production fiQ{ \T7\ 118] . Such calculations are not part of qcdnum itself, but can be 
coded in add-on packages; see Section [HI 

An important feature of VFNS evolution is that the input scale /ig cannot be above the 
lowest heavy flavour threshold /j,"^. This is because otherwise heavy flavour contributions 
must be included in the input parton densities which clearly is in conflict with the 
dynamic generation of heavy flavour by the QCD evolution equations. 
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Another feature of the VFNS is the existence of discontinuities at the flavour thresholds 
in as and the parton densities; we will now turn to the calculation of these discontinuities. 
Because the beta functions (12.21) depend on nf, it follows that the slope of the ag 
evolution is discontinuous when crossing a threshold in the VFNS. Beyond LO there are 
not only discontinuities in the slope but also in ctg itself. In N^LO, the value of Os"^ 
is, at a flavour threshold, related to Os"^ by [TT], [T^ 

ai"^^'\K^^l) = at^\K^il) + Y,{ y^'\^^l4.)Y Y.^n,\r.' «:} ^ = 1,2. (2.33) 

n=l j=0 

Here /i| is the threshold defined on the factorisation scale and k is the ratio f^'^/^p at fi\. 
For Gg = as/47r, the coefficients C in (I2.33P read 

/^ n r^ 2 ^ 14 ^ 38 r^ 4 

L'l.O — U, Oi^i — g, 02,0 — y, t^2,l — —, 02,2 — g- 

Note that there is always a discontinuity in as at NNLO. At NLO, a discontinuity 
only occurs when k 7^ 1, that is, when the renormalisation and factorisation scales are 
different. In case of upward evolution, as"^ is computed directly from (I2.33P while 
for downward evolution, as'^''^ is evaluated by numerically solving the equation 



where the function Aas(as) is given by the second term on the right-hand side of (I2.33p . 

In the VFNS at NNLO, not only as but also the parton densities have discontinuities at 
the flavour thresholds 



g{x,fil,nf + l) = g{x,fil,nf) + Ag{x,fil,nf) 

<lt{x,^J'l,nf + l) = qf{x,fil,nf) + Aqf{x,fil,nf) i = l,...,nf 

/i+(x,/i^,n/ + l) = Ah'^{x,iJ,l,nf) 

h-{x,fjl,nf + l) = Ah-{x,fxl,nf) =0, (2.34) 

where h = (c,b,t) for nj = (3,4,5). Note that a heavy quark h becomes a light quark 
qi above the threshold fi\. 

In QCDNUM, the flavour thresholds on the renormalisation scale are adjusted such that 
Uf changes by one unit in both the beta functions and the splitting functions when 
crossing a threshold. With this choice, the parton densities are continuous at LO and 
NLO while at NNLO the calculation of the discontinuities is considerably simplifled (all 
terms proportional to powers of ln(m^//i^) in ref. |20j vanish). So we may write 

^9{x,fJ'l,nf) = al{[Agq^qg]{x,fil,nf) + [Agg®g]{x,fil,nf)} 

A/i+(x,/i^,n/) = al{[Ah^®qs]{x,fil,nf) + [Ahg®g]{x,fil,nf)}. (2.35) 

Here a^ stands for as"'^ ('^/^a) ^^ deflned by (12.331) . The convolution kernels Aij can be 
found in Appendix B of [20] 



6ln the notation of HO], ^gq = ^q.?] (eq. B.5), Agg = A^f^ (B.7), Aq = A'^^;^^^ (B.4), Ai,q 
^Hq (B-1) and A^g — A^^ (B.3). For the latter we use a parametrisation provided by A. Vogt. 
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The discontinuities in the basis vectors {e^ ) are calculated from 
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ef (a;, fit n/ + 1) = ef (x, /i^, nj) + Aef{x, /i^, nj) + Xi{nf)Ah^{x, /i^, nj), (2.36) 

where the light component Ae^ is given by (12.351) . with q- replaced by e^ . With the 
definition (12.221) of the basis functions, the values of the coefficients Aj(nj-) are 



(2.37) 



When the densities are evolved upward in fi"^, it is straight forward to calculate with (I2.34p 
and (I2.35P the parton densities at nf + 1 from those at rif. However, qcdnum is capable 
to invert the relation between Uf and nj + 1 so that it can also calculate the disconti- 
nuities in case of downward evolution. For this it is convenient to write the calculation 
of the singlet and gluon discontinuities in matrix form, similar to (12. 9p 

In Section 13.31 we will show how (I2.38P is turned into an invertible matrix equation. 

Note that the heavy quark non-singlets do not obey the DGLAP evolution equations 
over the full range in /x^, because the heavy fiavours are simply set to zero below their 
thresholds, instead of being evolved. The evolution of the set \ef) thus proceeds in 
the VFNS as follows: The singlet /valence densities |e^) and the light non-singlets \ef^) 
are evolved both upward and downward starting from some scale /Iq < /i^. The heavy 
non-singlets je^gg) are dynamically generated from the DGLAP equations by upward 
evolution from the thresholds /i^bt- ^^ ^^^ below the thresholds, |e4 5g) is set equal 
to the singlet and |ejgg) to the valence. This is equivalent to setting the heavy quark 
and anti-quark distributions to zero, except that at NNLO the heavy fiavours do not 
evolve from zero but from the non-zero discontinuity given in (I2.34p . This is illustrated 
in Figure [2] where we plot the charm and bottom starting distributions, normalised to 
the singlet distribution. It is seen that the bottom discontinuity is less than 3% of the 
singlet over the whole range in x, while for charm it is much larger, exceeding 10% at 
low X. Note that the starting distributions are negative below x ~ 10~^. 



3 Numerical Method 

The DGLAP evolution equations are in qcdnum numerically solved on a discrete nxm 
grid in x and /i^. In such an approach the convolution integrals can be evaluated 
as weighted sums with weights calculated once and for all at program initialisation. 
Because of the convolutions, the total number of operations to solve a DGLAP equation 
is quadratic in n and linear in m. The accuracy of the solution depends, for a given 
grid, on the interpolation scheme chosen (linear or quadratic). 

The advantage of this 'x-space' approach, compared to W-space' [IH], is its conceptual 
simplicity and the fact that one is completely free to chose the functional form of the 
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Figure 2: The NNLO starting densities q^{x,^\), normalised to the singlet density qs{x,ij1), for 
charm (full curve) and bottom (dotted curve). 

input distribution since it is fed into the evolution as a discrete vector of input values. 
A disadvantage is that accuracy and speed depend on the choice of grid and that each 
evolution will yield no less than n x m parton density values (typically 10^) whether 
you want them or not. 

The numerical method used in qcdnum is based on polynomial spline interpolation 
of the parton densities on an equidistant logarithmic grid in x and a (not necessarily 
equidistant) logarithmic grid in /i^. The order of the x- interpolation can in be set to 
k = 2 (hnear) or 3 (quadratic). The interpolation in /i^ is always quadratic. With such 
an interpolation scheme, the DGLAP evolution equations transform into a triangular set 
of linear equations in the interpolation coefficients. This leads to a very fast evolution 
of these coefficients from some input scale /Xq to any other scale fif on the grid. In the 
following sections we will describe the spline interpolation, the calculation of convolution 
integrals and the QCD evolution algorithm. Note that several features of the qcdnumI? 
numerical method have been previously proposed in, for example, [211 [22] ■ 



3.1 Polynomial Spline Interpolation 

To interpolate a function /i(|/)|j we sample this function on an (n + l)-point grid 



yo<yi < ... < Vn-l < Vn 

and parametrise it in each interval by a piecewise polynomial of order k. Such a piecewise 
polynomial is turned into a spline by imposing one or more continuity relations at each 
of the grid points. Usually — but not always — continuity is imposed at the internal grid 
points on the function itself and on all but the highest derivative, which is allowed to be 



^In QCDNUM, h{y) represents a parton momentum density in the scaling variable y ^ —In a;. How- 
ever, for this section the identification of h with a parton density is not so relevant. 
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discontinuous. Without further constraints at the end points, the sphne has k + n — 1 
free parameters. Increasing the order k of the interpolation thus costs only one and 
not n extra parameters as is the case for unconstrained piecewise polynomials. 

It is convenient to write a spline function as a linear combination of so-called B-splines 



h{y) = J2Ay^iy)■ 



(3.1) 



The basis Yi of B-splines depends on the order k, on the distribution of the grid points 
along the y axis (equidistant in qcdnum) and on the number of continuity relations 
we wish to impose at the internal grid points and at the two end points. For how to 
construct a B-spline basis and for more details on splines in general we refer to [23j. 

In Figure [3] are shown the B-splines for linear {k = 2) quadratic {k = 3) and cubic 




Figure 3: B-spline bases generated on an equidistant grid, (a) Linear B-splines (fc = 2). Removing the 
dashed spline enforces the boundary condition h{yo) = 0; (b) Quadratic B-splines (fc = 3). Removing 
the first two dashed splines enforces the boundary condition h{yo) = h'{yo) = 0; (c) Cubic B-splines 
(fc = 4). Removing the first three dashed splines enforces the boundary condition h{yo) = h'{yQ) = 
h"{0) — 0. Spline interpolation on such a basis is numerically unstable. 

{k = 4) interpolation on an equidistant grid. In case /i(yo) = ^(0) = — which is always 
true for parton densities — we may remove the first B-spline in the plots of Figure |3l 
Removing the second B-spline in Figure |3]d gives quadratic interpolation with an addi- 
tional boundary condition h'(yo) = OJj With these boundary conditions — and because 



^A parton density parametrisation should thus behave like h{y — > 0) ex j/'^ with A > 1 because other- 
wise the condition h'{0) = is violated and the spline might oscillate. All known pdf parametrisations 
fulfil this requirement but when the parameters are under control of a fitting program one should take 
precautions that A will always stay above unity. 
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the grid is equidistant — the remaining B-sphnes possess translation invariance, that is, 
the basis can be generated by successively shifting the first spline one interval to the 
right (full curves in Figure |3^,b). Translation invariance greatly simplifies the evolution 
algorithm, as we will see later. 

It is therefore tempting to extend the scheme to cubic interpolation by removing the 
first three B-splines in Figure |3t. This would yield a translation invariant basis with 
the boundary conditions /i(yo) = h'{yo) = h"{yQ) = 0. However, it turns out that such 
a cubic spline interpolation tends to be numerically unstable. The cure is to drop the 
constraint /i"(z/o) = and impose a constraint on h'{yn) at the other end of the grid. But 
this does not fit in the evolution algorithm as it now stands so that we have abandoned 
cubic and higher order splines in qcdnum. 

If we number the B-splines 1,2, ... ,n from left to right as indicated in Figure [3] it is 
seen that for both k = 2 and 3 the following relation holds (translation invariance): 



F,(l/) = Fi(y-t/,_i). 
Furthermore, for linear interpolation {k = 2) we have Yi[yi) = 1 so that 



(3.2) 



h{yo) 







= AiY,{yi) = A, l<i<n. (3.3) 

Likewise, for quadratic interpolation (/c = 3) we have Yi_i{yi) = Yi{yi) = 1/2 so that 

hiyo) = 

%i) = A,Y,iy^) = A^/2 

h{yi) = A,^,Yi^iiyi) + A,Yiiyi) = {Ai_i + Ai)/2 2<%<n. (3.4) 

We denote h{%ji) by hi, the column vector of function values by /i = {h\, . . . , hn)^, the 
corresponding vector of spline coefficients by a and write (13. 3p and (13. 4p as 



h = Sa 



(3.5) 



where S is the identity matrix in case of linear interpolation and a lower diagonal band 
matrix for the quadratic spline. On a 5-point equidistant grid yo, . . . ,y4,, for instance, 
we have in case of quadratic interpolation the vector h = {hi, . . . , h^^j^ and the matrix 



/I 



\ 



V 



1 1 
1 1 

1 1/ 



with inverse 



1-1 



/ 1 

-1 
1 

V-1 



\ 



1/ 



(3.6) 



Note that S is sparse but S~^ is not. Thus, when a parton distribution Hq is given 
at some input scale /Iq, the corresponding vector ao of spline coefficients is found by 
solving (I3.5p rl This vector is then evolved to other values of jj? using the DGLAP 
evolution equations as is described in the next two sections. 



^Obtaining a from solving p.Sp by forward substitutfon (Appendix [BJ costs 0(27i) operations. This 
is cheaper than the alternative of calculating a = S^ h which costs 0(n^/2) operations. 
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3.2 Convolution Integrals 

The Mellin convolution fl2.5p calculated in qcdnum is not that of a number density / 
and some kernel g but, instead, that of a momentum density p = xf and a kernel q = xg. 
These convolutions differ by a factor x: 

[p^q]{x)=x[f^g]{x). (3.7) 

This also true for multiple convolution: for p = xf, q = xg and r = xh we have 

[p ^ q ^ r]{x) = x[f ® g ^ h]{x) . (3.8) 

A change of variable y = — In x turns a Mellin convolution into a Fourier convolution: 

ry rv 

[f ® g\{x) = [u ® "^Mv) = / (^z u{z) v{y — z) = I dz u{y — z) v{z), (3.9) 

Jo Jo 

where the functions u and v are defined by u{y) = /(e~^) and v{y) = g{e~'^). 

In the following we will denote by h{y, t) a parton momentum density in the logarithmic 
scaling variables y = —In a; and t = In/i^. In terms of h, the DGLAP non-singlet 
evolution equation (I2.12p is written as 

^^^ = J' dz Qiz, t) h{y -z,t) = j' dz Q{y - z, t) h{z, t) (3.10) 

with a kernel Q{y,t) = e^^P(e~^,t). Here P{x,t) is a non-singlet splitting function, as 
given in Section 12.21 To solve (I3.10p we first have to evaluate the Fourier convolution 

I{y,t)= [ dzQ{y-z,t)h{z,t). (3.11) 

Jo 

Inserting (13. ip in (13. lip we find for the integrals at the grid points yi (for clarity, we 
drop the argument t in the following) 



liVi) = V Aj r dz Q{y, - z) Yj{z) = V W.jAj (1 < i < n). 
.=1 ^0 



(3.12) 



The summation is over the first i terms only, because B-splines with an index j > i are 
zero in the integration domain z < yi, see Figure El 

Eq. (I3.12P defines the weights Wij which are calculated as follows. Because Yj{y) = 
for y < yj^i the weights can be written as 

Wij = dz Q{yi - z)Yj{z) = dz Q{yi - y,_i - z)Y^{z) (3.13) 

Jyj-i Jo 

where we have used (13. 2p in the second identity. From the property of equidistant grids 

Vi + Vj = Vi+j 
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it follows that Wij depends only on the difference i — j (Toeplitz matrix): 

Wij = Wi^j+i with Wf^= j dzQ{yt-z)Yi{z) (1 < ^ < n). (3.14) 

Jo 

The integrand only contributes in the region fcA where Yi is non-zero so that in practical 
calculations the upper integration limit yi is replaced by niin(?/£, /cA), with A the grid 
spacing. We remark that the calculation of the weights wi is a bit more complicated 
than suggested by (13.141) because singularities in the splitting functions have to be taken 
into account; for the relevant formula's we refer to Appendix \K\ 

The weights can thus be arranged in a lower-triangular Toeplitz matrix, as is illustrated 
by the 4x4 example below: 



W., 



W2 Wi 
W3 W2 Wi 
\ Wi W3 W2 Wi j 



(3.15) 



This matrix is fully specified by the first column, taking n instead of nin + l)/2 words 
of storage. This is not only advantageous in terms of memory usage but also in terms of 
computing speed since frequent calculations like summing the perturbative expansion 

W{t) = a,(t){ iy(°) + a,{t)W'-'^ + ■■■} (3.16) 

takes only 0(n) operations instead of 0(n^/2). We write the vector of convolution 
integrals as / and express 03.121) in vector notation as 

I=Wa. (3.17) 

Also multiple convolutions can be calculated as weighted sums. Let f{x) be a number 
density and Ka.b{x) be two convolution kernels. The vector of Mellin convolutions 

li = Xi[f Ka KbliXi) 

can be calculated from (I3.17p . using the weight table 

W = WaS-^Wb. (3.18) 

Here Wa and Wb are the weight tables of Ka and Kb, respectively, and S is the trans- 
formation matrix defined by (13. 5p . 

Another interesting convolution is that of two number densities fa and fb 

h = Xi[fa® fb\{Xi). 

This 'parton luminosity' [H] (times x) is calculated from the Fourier convolution 

I{y^)= / dzha{z)hb{yi-z). (3.19) 

Jo 

19 



Inserting the spline representation (13. ip gives an expression for the convolution integral 
as a weighted sum over the set of spline coefficients a of ha and b of hb, 

I I „y. 

^(Vi) = Y^Y^ ^J^k Wijk with Wijk = / dz Yj{z)Yk{yi - z). 
j=i k=i J^ 

To reduce the dimension of Wijk-, we use the translation invariance fl3.2p and write 

/•J/i-j + l 

Wijk= dzYi{z)Yk{yi-j+i- z). 

Jo 

Because B-splines with index k > z— j + 1 do not have their support inside the integration 
domain, we obtain an upper limit k < i—j + 1. Again using translation invariance yields 

ryi-j-k+2 
Wijk = dz Yi{z) Yi(yi_j_fc+2 - z). 

Jo 

We now have a compact expression for the convolution integral (J3TT9J): 

'' ^ ^zij^ rye 

HVi) = X] X^ ^j^k Wi_j^k+2 with we= dz Yi{z) Yi{yi - z). (3.20) 
j=i fc=i J^ 

Because Yi has a limited support, it turns out that only the first 3 (5) terms of w^ are 
non-zero in case of linear (quadratic) interpolation. The operation count to calculate a 
convolution of parton densities is thus not more than 0(5?2), for quadratic splines. 



3.3 DGLAP Evolution 

We denote by the vector Hq a non-singlet quark density at the input scale to = In /ip- The 
derivative of ho with respect to the scaling variable t is given by the DGLAP evolution 
equation (13.101) which can be written in vector notation as, from (13. 5p and (I3.17p 

— — = ——- = Wo ao or — — = a = S Wo ao- (3.21) 

dt dt dt 

Likewise we have at ti 

a[ = S-^Wi ai. (3.22) 

We have indexed the weight matrices above by a subscript because they depend on t 
through multiplication by powers of Og, see (I3.16p . 

Assuming that a{t) is quadratic in t, we can relate ao, o-i, (^o ^-nd a'^ by 

ai = ao + (a[, + a'JAi (3.23) 

with Ai = (ti — to)/2. If ti > to, Ai is positive and we perform forward evolution. If 
ti < to, Ai is negative and we perform backward evolution. 

Inserting (I3.2ip and (13.221) in (I3.23P we obtain a relation between the known spline 
coefficients ao and the unknown coefficients a^ 

(1 - 5"^WiAi) ai = (1 + 5"^WoAi) ao. (3.24) 
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Multiplying both sides from the left by Ui = S/Ai gives 

(C/i - W^) ai = (C7i + Wo) ao. (3.25) 

Eq. fl3.25p is more convenient than fl3.24p because matrix multiplication S'^W is re- 
placed by matrix addition]^ Note that C/ is a lower diagonal band matrix so that 
U ±W is still lower triangular with, in fact, the Toeplitz structure fl3.15p preserved. 
All this leads to a very simple and fast evolution algorithm, starting from ag: 

1. At to) calculate ao from fl3.5p . Wq from fl3.16p and t/i as defined above. Then 
construct the vector 61 = (C/i + Wo) ao- 

2. Subsequently, at ti, 

(a) Calculate Wi and the lower triangular matrix Vi = C/i — Wi; 

(b) Solve the equation Viai = bi by forward substitution, see Appendix [B| 

(c) Calculate C/2 and ^2 = (C/i + C/2)f*i — 61 for the next evolution to ^20 

3. Repeat step 2 at ^2 and so on. 

With this algorithm each evolution step consists of a few vector manipulations which 
have an operation count 0(n) and solving one triangular matrix equation which has an 
operation count 0(n^/2). The total operation count only very weakly depends on the 
order k of the interpolation chosen: quadratic interpolation is almost for free. 

The algorithm can also be used for the coupled evolution of the singlet quark (ttg) 
and gluon (ttg) spline coefficients, provided we make the following replacements in the 
formalism: 



aj \ SJ V^gq ^g 

In Appendix [B] is shown how the coupled triangular equations are solved by extending 
the forward substitution algorithm. The operation count is 4 x 0[n'^/2) so that for m 
grid points in t we have in total 0{2n'^m) operations for the singlet-gluon evolution and 
0{'n?m/2) operations for each non-singlet evolution. 

Finally, let us express in vector notation the NNLO parton density discontinuities at 
the flavour thresholds. The relation between the singlet and gluon distributions at Uf 
and n/ -(- 1 as given by f l2.38p can be written as 



^ J \"g/ V ^gq ^ + ^gg/ \"g/ 

It is easy to solve this linear equation for a("'/+^) when a*^"/) is known (forward evolution) 
or for a^"^-* when a*^"'-^"'"^-* is known (backward evolution). Likewise, we may write for 
the non-singlet discontinuities 

S at'^'^ = (5 + Aqq) at''^ + A (^hq at'^ + Ahg a^^) , (3.27) 



"'^^In fact, adding a matrix with band structure p.6p to a lower triangular matrix with structure (I3.15P 
takes only two additions irrespective of the dimension of the matrices. 

^^Using p.25p it is a simple exercise to establish this relation between b, U and a. Note that b in 
step (2c) is calculated much faster than b in step (1). 
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where A is defined by f l2.36p . Also this equation can easily be inverted. 

It can be seen from (I3.5p and fl3.23p that h{y,t) is, by construction, a spline in both 
the variables y and t. However, it turns out that it is technically more convenient to 
represent the pdfs by their values on the grid, instead of by their spline coefficients. 
Polynomial interpolation of order k in y and quadratic in t is then done locally on a 
k X 3 mesh around the interpolation point. The NNLO discontinuities are preserved 
by storing, at the flavour thresholds, the pdf values for both Uf — 1 and Uf, and by 
prohibiting the interpolation mesh to cross a flavour threshold. Note, however, that the 
interpolation routine yields a single-valued function of t, so that one has to calculate 
h{y, tc,b,t — e) to view the discontinuityo 

In QCDNUM it is possible to evolve on multiple equidistant y-grids which allow for a 
flner binning at low y (large x) where the parton densities are rapidly varying. This is 
illustrated below by a grid Gq which is built-up from three equidistant sub-grids Gi, G2 
and G3 with spacing A/4, A/2 and A, respectively. 



Go 

Gi 
G2 
G3 



Z/o 



yi 



1/2 



2/3 



(I) 


(11) 


(III) 






A 



On such a multiple grid, the parton densities are flrst evolved on the grid Gi and the 
results are copied to the region (I) of Go- The evolution is then repeated on the grids G2 
and G3 followed by a copy to the regions (II) and (III) of Go, respectively. We refer to 
Section 14.31 for spectacular gains in accuracy that can be achieved by employing these 
multiple grids. 

As remarked above, the evolution algorithm can — at least in principle — handle both 
forward and backward evolution in /i^ simply by changing the sign of A in (I3.23p . This 
works very well for linear spline interpolation but it turns out that backward evolution 
of quadratic splines can sometimes lead to severe oscillations. This is illustrated in 
Figure H] where is shown a non-singlet quark density evolved downward from /Xq = 5 
to /x^ = 2 GeV^ in the quadratic interpolation scheme (dotted curve). In qcdnum this 
numerical instability is handled as follows: (i) evolve downward from Hq to n^ in the 
linear interpolation scheme (which is stable); (ii) then take /i^ as the starting scale and 
evolve upward to /ig in the quadratic interpolation scheme (also stable); (iii) calculate 
the difference A/ between the newly evolved pdf and the original one at /ig; (iv) subtract 
A/ from the starting value at /ig used in (i) and repeat the procedure. 

The full curve in the top plot of Figure H] shows the result of downward evolution in the 
linear interpolation scheme, that is, without iterations. Oscillations are absent but the 
evolution is not very accurate as is evident from the difference between the dotted and 
full curve at large x. One iteration already much improves the precision as can be seen 
from the good match at large x between the two curves in the bottom plot. It turns out 
that one iteration (qcdnum default), perhaps two, are sufficient while more iterations 



Do not take e too small because qcdnum may snap to the threshold, see Section 15^2 
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Figure 4: A non-singlct parton density xf(x) versus log(x) evolved downward from /Xq = 5 to /i^ = 
2 GeV^ in the quadratic interpolation scheme showing large oscillations (dotted curve). The full curve 
in the top plot shows the result of downward evolution in the linear interpolation scheme. The full 
curve in the bottom plot shows an improved result obtained by iteration, as described in the text. 



tend to spoil the convergence. Clearly best is to limit the range of downward evolution 
by keeping /Iq low or, if possible, to set it at the lowest grid point to avoid downward 
evolution altogether. 

QcDNUM checks for quadratic spline oscillation as follows. We denote the values of 
the quadratic B-spline at (|A, A, |A) by (61,62,^3) = (|, \, f)- It is easy to show that 
quadratic interpolation mid-between the grid points is given hy u = Da, where D is 
a lower diagonal Toeplitz band matrix, of bandwidth 3, which is characterised by the 
vector (61, 63, hi). Likewise, the linear interpolation of the spline at the mid-points is cal- 
culated from V = Ea, where E is the lower diagonal Toeplitz band matrix (I&25 &2, 2^2)- 
The maximum deviation e = \\u — v\\ = \\{D — E)a\\ should be small; for pdfs sampled 
on a reasonably dense grid, e ~ 0.1 or less. For each pdf evolution, e is computed at the 
input scale, and at the lower and upper end of the /x^ grid. An error condition is raised 
when it exceeds a given limit, indicating that the spline oscillates, or that the x-grid is 
not dense enough. 
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4 The QCDNUM Program 

4.1 Source Code 

The QCDNUM source code can be downloaded from the web site 

http : //www . nikhef . rLl/user/h24/qcdnuni 

Unpacking the tar file produces a directory qdcnum-xx-yy with xx-yy the version num- 
ber. Sub-directories contain the source code, example jobs, write-up and a simple script 
to make a QCDNUM library, see the README file. The code comes with a utility package 
MBUTIL (including write-up) which is a collection of general-purpose routines (some de- 
veloped privately, some taken from CERNLIB and some taken from public source code 
repositories like netlib). Because QCDNUM uses several of these routines, mbutil must 
also be compiled and linked to your application program. Apart from this, QCDNUM is 
completely stand-alone. To calculate structure functions, the zmstf and hqstf add-on 
packages are provided, see Sections IC.3I and ID.li 

Before compiling QCDNUM you may want to set several parameters which control the 
size of internal arrays. These parameters can be found in the include file qcdnum. inc: 

mxgO Maximum number of multiple x-grids [5] . 

mxxO Maximum number of points in the x-grid [300]. 

mqqO Maximum number of points in the /i^-grid [150]. 

mptO Maximum number of interpolations calculated in a single call [5000]. 

miwO Maximum number of information words in a weight store [20] . 

mbf Maximum number of fast convolution scratch buffers [20]. 

nwf Size of the qcdnum dynamic store in words [400000]. 

The first 6 parameters are simply dimensions of book-keeping arrays which you may 
want to adjust to your needs. More important is the parameter nwf that defines the 
size of an internal store that contains the weight tables and the tables of parton densities. 
How many words are needed depends on the size of the tables which, in turn, depends 
on the size of the current x-n"^ grid. It also depends on how many different sets of tables 
(un-polarised pdfs, polarised pdfs, fragmentation functions, etc.) one wants to store. In 
this respect, qcdnum is very user-friendly by always gracefully grinding to a halt if it 
runs out of memory, with a message that tells how large nwf should be. 

4.2 Application Program 

To illustrate the use of QCDNUM, we present in FigureOthe listing of a simple application 
program. For a detailed description of the subroutine calls, and for additional routines 
not included in the example, we refer to Section [51 

The first step in a qcdnum based analysis is initialisation (qcinit), setting up the a;-/i^ 
grid (gxmake, gqmake) and the calculation of the weight tables (f illwt). The weights 
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c 
c 




program example 








impl; 


Lcit double precision (a 


-h,o-z) 










data 


ityp/1/, iord/3/, nfin/0/ 




lunpolarised, NNLO, VFNS 






data 


asO/0.364/, r20/2.D0/ 






! alphas 






external func 






! input part on dists 






dimension def (-6:6,12) 






! flavor decomposition 






data 


def / 








c— 




tb bb cb sb ub db g 


d u 


s 


c b t 


c— 




-6 - 


-5 -4 -3 -2 -1 


1 2 


3 


4 5 6 




+ 


0., 


0. , 0., 0., 0.,-l., 0. , 


1. , 0. , 


0., 


. , . , . , ! dval 




+ 


0., 


0. , 0., 0.,-l., 0. , 0. , 


0. , 1. , 


0., 


. , . , . , ! uval 




+ 


0., 


0. , 0.,-l., 0., 0. , 0. , 


0. , 0. , 


1., 


. , . , . , ! sval 




+ 


0., 


0. , 0., 0., 0., 1. , 0. , 


0. , 0. , 


0., 


. , . , . , ! dbar 




+ 


0., 


0. , 0., 0., 1., 0., 0. , 


0., 0., 


0., 


. , . , . , ! ubar 




+ 


0. , 


0. , 0., 1., 0., 0. , 0. , 


0., 0., 


0., 


. , . , . , ! sbar 




+ 


78*0. / 












data 


xmin/l.D-4/, nxin/100/, 


iosp/3/ 




!x grid, splord 






dimension qq(2),wt(2) 






!mu2 grid 






data 


qq/2.D0,l.D4/, wt/l.DO, 


l.DO/, nqin/60/ !mu2 grid 






data 


q2c/3.D0/, q2b/25.D0/, 


qO/2.0/ 




! thresh and mu20 


c 




data 


x/l.D-3/, q/l.D3/, qmz2/8315 . 25D0/ 


! output scales 


call 


qcinit (6, ' ' ) 






! initialise 






call 


gxmake (xmin , 1 , 1 , nxin , nx 


,iosp) 




!x-grid 






call 


gqmake ( qq , wt , 2 , nqin , nq) 






!mu2-grid 






call 


f illwt (ityp , idl , id2 ,nw) 






! calculate weights 






call 


setord(iord) 






!L0, NLD, NNLO 






call 


setalf (as0,r20) 






! input alphas 






iqc 


= iqfrmq(q2c) 






!mu2c 






iqb 


= iqfrmq(q2b) 






!mu2b 






call 


setcbt (nf in, iqc , iqb,0) 






! thresholds in the VFNS 






iqO 


= iqfrmq(qO) 






! starting scale 






call 


evolfgC ityp, func ,def , iqO,eps) 




! evolve all pdfs 






csea 


= 2.D0*fvalxq(ityp,-4,x 


,q,0) 




! charm sea at x,Q2 






asmz 


= asfunc(qmz2,nf out,ierr) 




! alphas (mz2) 


c 
c 




end 












double precision function func(id,x) 




! momentum density xf (x) 



implicit double precision (a-h,o-z) 
if (id.eq.O) func = gluon(x) 
if(id.eq.l) func = dvalence(x) 



! = always gluon 
!1 = defined in def 



if(id.eq.6) func = strangebar (x) 

return 

end 



!6 = defined in def 



Figure 5: Listing of a QCDNUM application program evolving a complete set of parton densities in the 
VFNS at NNLO. The array def defines the light quark valence (xq—xq) and anti-quark (xq) distributions 
as an input to the evolution. The x dependence of the input densities is coded in the function func. 
After evolution, the pdfs are interpolated to some x and ^^ and as{m'^) is calculated. 
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depend on the grid definition and the interpolation order so that f illwt must be called 
after the grid has been defined. The weight tables are calculated for LO, NLO and 
NNLO as well as for all possible flavour settings in the range 3 < n/ < 6 so that you 
do not have to call f illwt again when you set or re-set QCDNUM parameters further 
downstream. Although the weight calculation is fast (typically about 10-20 s) it may 
become a nuisance in semi-interactive use of qcdnum so that there is a possibility to 
dump the weights to disk and read them back in the next QCDNUM run. 

In the example code, the weight calculation is followed by setting the perturbative 
order (setord) and the input value of ttg at some renormalisation scale fi^ (setalf ). 
The call to setcbt sets the VFNS mode and defines the thresholds on the factorisation 
scale /ip. All the calls that set evolution parameters are destructive in the sense that 
they invalidate the parton densities in memory, if any. In this way all qcdnum results 
are consistently obtained with the same value of ag, the same perturbative order, etc. 

The second step is to evolve the parton densities from input specified at the scale /Iq. It 
is important to note that QCDNUM evolves parton momentum densities xf{x), although 
all theory in this write-up is expressed in terms of parton number densities f{x). The 
evolution is done by calling the routine evolfg which evolves 2n/ -|- 1 input parton 
densities (quarks plus gluon) in the FFNS or VFNS scheme. The routine internally takes 
care of the proper decomposition of the input quark densities into singlet and non- 
singlets. In the VFNS the input scale /Xq must lie below the charm threshold /i^ so that, 
as a consequence, /x^ must lie above the lower boundary of the /i^ grid. 

The flavour composition of each of the input quark densities is given by a table of 
weights def (-6:6, 12). In the example program, six light quark input densities are 
defined: three valence densities x{q — q) and three anti-quark densities xq. This is 
sufficient input to run evolutions in the VFNS scheme. One is completely free to define 
the flavour composition of the input quark densities as long as they form a linearly 
independent set (qcdnum checks this). Note that the flavours are ordered according to 
the PDG convention d, u, s, . . . and not u, d, s, . . . as often is the case in other programs. 

The X dependence of these momentum densities at /Iq must be coded for each identifier 
in an if-then-else block in the function func. The sum rules 



xg{x)dx + / 
Jo 


xqs{x)(\x 


[ [d{x) - 
Jo 


- d{x)](\x 


j\u{x)- 
Jo 


- u{x)\(\x 



: 2 (4.1) 

cannot be reliably evaluated by qcdnum since it has no information on the x-dependence 
of the pdfs below the lowest grid point in x. These sum rules should therefore be built 
into the parametrisation of the input densities. The evolution does, of course, conserve 
the sum rules once they are imposed at /Xq- The easiest way to evolve with a symmetric 
strange sea is to include xs^ = x{s — s) in the collection of input densities and set it to 
zero for all x at the input scale /Xq. In the VFNS at LO or NLO, the generated heavy 
flavour densities h = (c, b, t) are always symmetric {xh — xh = 0) but this is not true 
anymore at NNLO, which generates a small asymmetry. 

26 



After the parton densities are evolved, the results can be accessed by f valxq. This 
routine transforms the parton densities from the internal singlet /non-singlet basis to 
the flavour basis and returns the gluon, a quark, or an anti-quark momentum density, 
interpolated to x and /x^. Also here the flavours d, u, s, . . . are indexed according to the 
PDG convention. The last call in the example program evolves the input value of a^ to 
the scale m%. This evolution is completely stand-alone and does not make use of the /i^ 
grid. The function asf unc can thus be called at any point after the call to qcinit. We 
refer to Section [5] for more ways to access the qcdnum results, and for ways to change 
the renormalisation scale with respect to the factorisation scale. 

Qcdnum has an extensive checking mechanism which maintains internal consistency 
and verifies that all subroutine arguments supplied by the user are within their allowed 
ranges. Error messages might pop-up unexpectedly when the renormalisation scale is 
changed with respect to the factorisation scale because the low end of the /i^ grid may 
then map onto values of /ip^ < A^ . 

Another QCDNUM feature is that nj = (4,5,6) and not (3,4,5) at the heavy flavour 
thresholds yU^. This implies, first of all, that parton evolution in the VFNS must start 
from /Xq < fj^l and not from /Xq < /x^, simply because the number of flavours must be 
Hf = 3 at the starting scale. There is, however, no restriction on the starting (renor- 
malisation) scale of as so that it may very well coincide with a flavour threshold, either 
before or after varying the renormalisation scale with respect to the factorisation scale. 
If this happens at NNLO, the input value of Os is assumed to include the discontinuity. 



4.3 Validation and Performance 

The CPU time that is needed to evolve a pdf on a discrete grid grows quadratically 
with the number of grid points in x. With linear (quadratic) interpolation the accuracy 
increases linearly (quadratically) with the number of grid points. It follows that an 
r-fold gain in accuracy will cost a factor of r^ in CPU for linear interpolation but only 
a factor of r for quadratic interpolation. This reduction in cost motivated the inclusion 
of quadratic splines in qcdnum. 

To investigate the performance of the two interpolation schemes, we compare results 
from QCDNUM to those from the A^-space evolution program PEGASUS [TU]. In this 
comparison a default set of initial distributions [25] is evolved at NNLO from /i^ = 2 
to yU^ = 10^ GeV^ with Uf = A flavours. The dashed curve in the top plot of Figure [6] 
shows the relative difference Ag/g versus x for qcdnum evolution with linear splines 
on a single 200 point grid extending down to s = 10~^. The accuracy at low x is 
satisfactory (few permille) but deteriorates rapidly to Ag/g > 2% for x > 0.35. 

The precision is much improved by evolving on multiple grids (Section 13. 3p as shown by 
the full curve in the top plot of Figure O Here the 200 grid points are re-distributed over 
flve sub-grids with lower limits as given in Table [TJ For each successive grid the point 
density is twice that of the previous grid. It is seen from Figure [6] that the precision is 
now better than 2% for x < 0.85. 

The dotted curves in Figure [6] (top and bottom) correspond to evolution with quadratic 
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0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 X 



Figure 6: The relative difference /S.g/g (in percent) of gluon densities evolved from /i^ = 2 to /i^ = 
10'* GeV^ by QCDNUM and PEGASUS. Top: Evolution with Uncar splines on a 200 point single grid 
down to a: = 10^^ (dashed curve) and on multiple grids (full curve). Bottom: Evolution with quadratic 
splines on a 100 point single grid (dotted curve, also shown in the top plot) and on multiple grids (full 
curve). Note the different vertical scales in the two plots. 



splines on a single 100 point grid. There is a large improvement in accuracy (more 
than a factor of 10) compared to linear splines even though the number of grid points is 
reduced from 200 to 100. However, also here the precision deteriorates with increasing x, 
reaching a level of 2% at x = 0.65. A five-fold multiple grid with lower limits as 
listed in Table [T] yields a precision Ag/ g < 5 x 10~^ over the entire range x < 0.9 as 
can be seen from the full curve in the lower plot of Figure O Note that this is for 
evolution up to /i^ = 10"^ GeV^; at lower /i^ the accuracy is even better since it increases 
(roughly linearly) with decreasing ln(/i^). To fully validate the qcdnum evolution with 
PEGASUS^fl we have made additional comparisons in the FFNS with nj = 3, 5 or 6 
flavours, in the VFNS with and without backward evolution, and with the renormalisation 
scale set different from the factorisation scale. This for both un-polarised evolution up 
to NNLO and polarised evolution up to NLO. 

As remarked in Section [X^ the quadratic spline evolution is not more expensive in CPU 
time than linear spline evolution. On the contrary: QCDNUM runs 4 times faster since 



-'^^Similar benchmarking between hoppet [21 and PEGASUS is given in pS] and |551, where also pdf 
reference tables can be found. We do not provide here benchmark tables for qcdnum, but a program 
that generates such tables and compares them with PEGASUS is available upon request from the author. 
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Table 1: Lower x limits of multiple grids used in the evolution with linear and quadratic splines. 

n Xi X2 X3 X4 x^ 



Linear interpolation 


200 


10-5 


0.01 


0.10 


0.40 


0.70 


Quadratic interpolation 


100 


10-5 


0.20 


0.40 


0.60 


0.75 


Relative point density 




1 


2 


4 


8 


16 



we need only 100 instead of 200 grid points. With the multiple grid definition given in 
Table [1] for quadratic sphnes, the density of the first grid (x > 10-^) is 12 points per 
decade. It follows that for evolution down to x = 10"^ (10"^) a grid with 100 + 12 = 112 
(100 — 12 = 88) points should be sufficient. 

To investigate the execution speed we did mimic a QCD fit by performing 1000 NNLO 
evolutions in the VFNS (13 pdfs), using a 60 point /i^ grid and the 5-fold 100 point x-grid 
given in Table [TJ After each evolution, the proton structure functions F2 and Fl were 
computed at NNLO for 1000 interpolation points in the HERA kinematic range. For this 
test, QCDNUM, MBUTIL, and ZMSTF were compiled with the gfortran compiler, using 
level 2 optimisation and without array boundary check. The computations took 18.5 
CPU seconds on a 2 GHz Intel Core 2 Duo processor under Mac OS-X: 8.5 s for the 
evolutions and 10 s for the structure functions. 



5 Subroutine Calls 

In this section we describe all available qcdnum subroutines and functions. For conve- 
nience a list of these is given in Table [H In the following we will prefix output variables 
with an asterisk (*). We use the FORTRAN convention that integer variable and function 
names start with the letters I-N. Character variables are given in quotes as in ' opt ' . 
Other variables and functions are in double precision unless otherwise stated. Note that 
fioating point numbers should be entered in double precision format: 



ix = ixfrmx ( x ) 

ix = ixfrmx ( O.IDO ) 

ix = ixfrmx (0.1 ) 



ok 
ok 
wrong ! 



Unlike FORTRAN, qcdnum is case insensitive so that character arguments like 'ALIM' 
or 'Alim' are both valid inputs. 

Most QCDNUM functions will, upon error, generate an error message. The inclusion of 
function calls in print or write statements can then cause program hang-up in case 
the function tries to issue a message. Thus: 

write(6,*) 'Glue = ', fvalxq(l,0,x,q, 1) ! not recommended 

glue = fvalxq(l,0,x,q,l) ! OK 

write (6,*) 'Glue = ', glue ! OK 
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Table 2: Subroutine and function calls in qcdnum. 



Subroutine or 


function 


Description 


QCINIT ( 


; lun, 'filename' ) 


Initialise 


SETLUN ( 


[ lun, 'filename' ) 


Redirect output 


SETIGETVAL ( 


; ' opt ' , val ) 


Set Get parameters 


SETIGETINT ( 


[ ' opt ' , ival ) 


Set Get parameters 


GXMAKE ( 


[ xmi, iwt, n, nxin, *nxout, lord ) 


Define x grid 


IXFRMX ( 


: X ) 


Get ix from x 


XFRMIX ( 


: ix ) 


Get X from i^ 


XXATIX ( 


; X, ix ) 


Verify grid point 


GQMAKE ( 


; qarr, wt, n, nqin, *nqout ) 


Define /i| grid 


IQFRMQ ( 


: q2 ) 


Get i^ from /ip 


QFRMIQ ( 


: iq ) 


Get /ip from i^ 


QQATIQ ( 


: q2, iq ) 


Verify grid point 


GRPARS ( 


; *nx, *xl, *x2, *nq, *ql, *q2, *io ) 


Get grid definitions 


GXCOPY ( 


[ *array, n, *nx ) 


Copy X grid 


GQCOPY ( 


[ *array, n, *nq ) 


Copy /i^ grid 


FILLWT ( 


[ itype, *idmi, *idma, *nw ) 


Fill weight tables 


FILLWC ( 


[ mysub, *idmi, *idma, *nw ) 


Custom weights 


DMPWGT ( 


; itype, lun, 'filename' ) 


Dump weight tables 


READWT ( 


; lun, 'fn', *idmi, *idma, *nw, *ie ) 


Read weight tables 


NWUSED ( 


; *nwtot, *nwuse, *nwtab ) 


Memory words used 


SETIGETORD ( 


; lord ) 


Set Get order 


SETIGETALF ( 


: alfs, r2 ) 


Set Get tts start value 


SETCBT ( 


' nfix, iqc, iqb, iqt ) 


Set Hf or thresholds 


GETCBT ( 


; *nfix, *q2c, *q2b, *q2t ) 


Get Hf or thresholds 


SETIGETABR ( 


; ar, br ) 


Set Get ;U^ scale 


RFROMF ( 


[ fscale ) 


Convert /i| to fi^ 


FFROMR ( 


' rscale ) 


Convert fi^ to /i| 


ASFUNC ( 


[ r2, *nf, *ierr ) 


Evolve as{n^) 


EVOLFG ( 


; itype, func, def, iqO, *eps ) 


Evolve all pdfs 


PDFINP ( 


; subr, iset, offset, *epsi, *nwds ) 


Pdfs from outside 


CHKPDF ( 


; iset ) 


True if pdf set exists 


FVALXQ ( 


; iset, id, X, qmu2, ichk ) 


Interpolate g, q, q) 


FVALIJ ( 


; iset, id, ix, iq, ichk ) 


\g,q,q) at grid point 


FPDFXQ ( 


[ iset, X, qmu2, *pdfs, ichk ) 


All pdfs 1^, g, q) 


FPDFIJ ( 


[ iset, ix, iq, *pdfs, ichk ) 


All pdfs \g, q, q) 


FSUMXQ ( 


[ iset, def, x, qmu2, ichk ) 


Linear combination 


FSUMIJ ( 


[ iset, def, ix, iq, ichk ) 


Linear combination 


FSNSXQ ( 


[ iset, id, X, qmu2, ichk ) 


Interpolate g, e^) 


FSNSIJ ( 


; iset, id, ix, iq, ichk ) 


\g, e^) at grid point 


FSPLNE ( 


; iset, id, X, iq ) 


Spline interpolation 


SPLCHK ( 


; iset, id, iq ) 


Check spline 



Output arguments are pre-fixed with an asterisk (*). 
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5.1 Initialisation 



call QCINIT ( lun, 'filename' ) 



Initialise qcdnum and define the output stream. Should be called before anything else. 

lun Output logical unit number. When set to 6, qcdnum messages appear 

on the standard output. When set to -6, the qcdnum banner printout 
is suppressed on the standard output. 

'filename' Output file name. Irrelevant when lun is set to 6 or -6. 



call SETLUN ( lun, 'filename' ) 



Redirect the qcdnum messages. The parameters are as for qcinit above. This routine 
can be called at any time after qcinit. 



call SETVAL I GETVAL ( 'opt', val ) 



Set or get qcdnum floating point parameters. 

'null' Result of a calculation that cannot be performed. Default, null = l.Dll. 
' epsi ' The tolerance level in the floating point comparison |x— ?/| < e, which qcdnum 

uses to decide if x and y are equal. Default, epsi = 1 .D-9. 
'epsg' Required numerical accuracy of the Gauss integration in the calculation of 

weight tables. Default, epsg = l.D-7. 
'elim' Allowed difference between a quadratic and a linear spline interpolation mid- 

between the grid points in x. Default, elim = 0.5; larger values may indicate 

spline oscillation. To disable the check, set elim < 0. 
'alim' Maximum allowed value of q;s(/^^)- When a^ exceeds the limit, a fatal error 

condition is raised. Default, alim = lOo 
'qmin' Smallest possible lower boundary of the /^^ grid. Default, qmin =0.1 GeV^. 
'qmax' Largest possible upper boundary of the /x^ grid. Default, qmax = 1 .Dll GeV^. 

These parameters can be set and re-set at any time after qcinit. 



call SETINT I GETINT ( 'opt', ival ) 



Set or get qcdnum integer parameters. 

' iter ' Set the number of iterations in the backward evolution. When set negative, one 
will evolve backward in the same interpolation scheme as the forward evolution 
(not recommended). When set to zero, one will evolve backward in the linear 
interpolation scheme, without iterations (not recommended either). A value 
larger than zero gives the number of iterations to perform. Default, iter = 1. 
This parameter is only relevant when one works with quadratic splines. 



-'^^When you raise alim > 10 then ag will at some point be limited by internal cuts in qcdnum. 
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'lunq' Retrieve the qcdnum logical unit number. Useful if one wants to write mes- 
sages on the same output stream as qcdnum. This option is only available for 
getint and not for set int. 

'ntab' Number of scratch buffers (maximum 20) generated by fastini (fast convolu- 
tion engine, Section 1^75]) . Will have no effect when set after the call to fastini. 
Default, ntab = 5. If one wants to generate more than 20 buffers, the value of 
mbfO in qcdnum . inc should be increased, and QCDNUM re-compiled. 

5.2 Grid Definition 

A proper definition of the grid in x and /i^ is important because it determines the speed 
and accuracy of the qcdnum calculations. The grid definition also governs the partition 
of the internal store which contains the weight tables and tables of parton densities. In 
addition, the routines set-up the bases of B-splines. 

The X grid must be strictly equidistant in the variable y = — In a; but in qcdnum one 
can generate multiple equidistant grids (Section 13. 3p to obtain a finer binning at low y 
(large x). Multiple grids are generated when the x-range is subdivided into regions with 
different densities, as is described below. 

The fi^ grid does not need to be equidistant. So one can either enter a fully user-defined 
grid or let qcdnum generate one by an equidistant logarithmic fill-in of a given set of 
intervals in /i^. 



call GXMAKE ( xmin, iwt, n, nxin, *nxout, lord ) 



Generate a logarithmic x-grid. 

xmin Input array containing n values of x in ascending order: xmin(l) defines the 
lower end of the grid while the other values define the approximate positions 
where the point density will change according to the values set in iwt. The list 
may or may not contain x = 1 which is ignored anyway. 

iwt Input integer weights. The point density between xmin(l) and xmin(2) will be 
proportional to iwt(l), that of the next region will be proportional to iwt (2) 
and so on. The weights should be given in ascending order and must always 
be an integer multiple of the previous weight. Thus, to give an example, the 
triplets {1,1,1} and {1,2,4} are allowed but {1,2,3} is not. 

n The number of values specified in xmin and iwt. This is also the number of 

sub-grids used internally by qcdnum. 

nxin Requested number of grid points (not including the point x = 1). Should of 
course be considerably larger than n for an x-grid to make sense. 

nxout Number of generated grid points. This may differ slightly from nxin because of 
the integer arithmetic used to generate the grid. 

lord One should set lord = 2 (3) for linear (quadratic) spline interpolation. 

With this routine, one can define a (logarithmic) grid in x with higher point densities 
at large x, where the parton distributions are strongly varying. Thus 
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xmin = l.D-4 

iwt = 1 

call gxmake (xmin , iwt , 1, 100, nxout, lord) 

generates a logarithmic grid with exactly 100 points in the range 10~^ < x < 1, while 

xmin(l) = l.D-4 

iwt(l) = 1 

xmin(2) = 0.7D0 

iwt (2) = 2 

call gxmake (xmin , iwt , 2 , 100, nxout, lord) 

generates a 100-point grid with twice the point density above x ~ 0.7. 
A call to gxmake invalidates the weight tables and the pdf store. 



ix = IXFRMX ( X ) X = XFRMIX ( ix ) L = XXATIX ( x, ix ) 



The function ixf rmx returns the index of the closest grid point at or below x. Returns 
zero if x is out of range (note that x = 1 is outside the range) or if the grid is not 
defined. The inverse function is x = xfrmix(ix). Also this function returns zero if ix 
is out of range or if the grid is not defined. To verify that x coincides with a grid point, 
use the logical function xxatix, as in 

logical xxatix 

ix = xfrmix(x) !x is at or above grid point ix 

if (xxatix(x,ix)) then !x is at grid point ix 

Note that qcdnum snaps to the grid, that is, x is considered to be at a grid point i if 
\y — yi\ < e with y = —Inx and, by default, e = 10~^. 



call GQMAKE ( qarr, wgt, n, nqin, *nqout ) 



Generate a logarithmic /ip grid on which the parton densities are evolvedo 

qarr Input array containing n values of /i^ in ascending order: qarr(l) and qarr(n) 
define the lower and upper end of the grid, respectively. The lower end of the 
grid should be above 0.1 GeV^. If n > 2 then the additional points specified in 
qarr are put into the grid. In this way, one can incorporate a set of starting 
values /Iq, or thresholds /^cbt- 

wgt Input array giving the relative grid point density in each region defined by qarr. 
The weights are not restricted by integer multiples as in gxmake but can be set 
to any value in the range 0.1 < w < 10. With these weights, one can generate 
a grid with higher density at low values of /i^ where as is changing rapidly. 
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Note that as is evolved (without using a grid) on /i^ which may or may not be different from /ip. 
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n The number of values specified in qarr and wgt (n > 2). 

nqin Requested number of grid points. If nqin < n then the grid is not generated 
but taken from qarr. This feature allows you to read-in your own /i^ grid. 

nqout Number of generated grid points. This may differ slightly from nqin because of 
the integer arithmetic used to generate the grid. 

A call to gqmake invalidates the weight tables and the pdf store. 



iq = IQFRMQ ( q2 ) q2 = QFRMIQ ( iq ) L = QQATIQ ( q2, iq ) 



The function iqf rmq returns the index of the closest grid point at or below /i^. The 
inverse function is qf rmiq. To verify that /i^ coincides with a grid point, use the logical 
function qqatiq. As described above for the corresponding x grid routines, a value of 
zero is returned if q2 and/or iq are not within the range of the current grid, or if the 
grid is not defined. 



call GRPARS ( *nx, *xnii, *xma, *nq, *qnii, *qma, *iord ) 



Returns the current grid definitions 

nx Number of points in the x grid not including x = 1. 

xmi Lower boundary of the x grid. 

xma Upper boundary of the x grid. Is always set to xma = 1. 

nq Number of points in the /i^ grid. 

qmi Lower boundary of the /i^ grid. 

qma Upper boundary of the fj,"^ grid. 

lord Order of the spline interpolation (2 = linear, 3 = quadratic). 



call GXCOPY ( *array, n, *nx ) 



Copy the x grid to a local array 

array Local array containing on exit the x grid but not the value x = 1. 

n Dimension of array as declared in the calling routine. 

nx Number of grid points copied to the local array. A fatal error occurs if array is 
not large enough to contain the current grid. 



call GQCOPY ( *array, n, *nq ) 



As above, but now for the /i^ grid. 

34 



5.3 Weights 

In this section we describe routines to calculate the weight tables, to dump these to disk 
and to read them back. The weight tables are calculated for all orders (LO,NLO,NNLO) 
and all number of flavours nj = (3, 4, 5, 6), irrespective of the current qcdnum settings. 
Tables can be created for un-polarised pdfs, polarised pdfs and fragmentation functions. 
All these pdf types can exist simultaneously in memory. For each type, one gluon table 
and 12 quark tables are generated by the routines, in addition to the weight tables. 



call FILLWT ( itype, *idniin, *idmax, *nwds ) 



Partition the pdf store and fill the weight tables used in the calculation of the convolution 
integrals. Both the x and fi"^ grid must have been defined before the call to f illwt. 

itype Select un-polarised pdfs (1), polarised pdfs (2) or fragmentation functions (3). 

Any other input value will select un-polarised pdfs (default), 
idmin Returns, on exit, the identifier of the first pdf table. Always idmin = 0. 
idmax Identifier of the last pdf table in the store. Always idmax = 12. 
nwds Total number of words used in memory. 

One can create more than one set of tables tables by calling f illwt with different values 
of itype. For instance, the sequence 

call f illwt (1, idmin, idmax, nw) lUnpolarised pdfs 
call f illwt (2, idmin, idmax, nw) ! Polarised pdfs 

makes both the un-polarised and the polarised pdfs available. For each pdf type, f illwt 
creates 13 pdf tables. If there is not enough space in memory to hold all the tables, 
f illwt returns with an error message telling how much memory it needs. One should 
then increase value of nwfO in the include file qcdnum. inc, and recompile QCDNUM '^^' 
Note that f illwt acts as a do-nothing when the pdf type already exists in memory: 

call f illwt (1, idmin, idmax, nw) lUnpolarised pdfs 
call f illwt (1, idmin, idmax, nw) !Do nothing 



call DMPWGK itype, lun, 'filename' ) 



Dump the weight tables (not the pdf tables) of a given pdf type to disk. When itype = 0, 
all pdf types in memory are dumpedo Fatal error if itype does not exist. Additional 
information about the qcdnum version, grid definition and partition parameters is also 
dumped, to protect against corruption of the dynamic store when the weights are read 
back in future QCDNUM runs. The dump is unformatted so that the output file cannot 
be exchanged across machines. 



^^This one may have to repeat several times, because f illwt proceeds in stages and is ignorant of 
the memory requirements of the next stage. 

^"^ This does not include the weight tables of custom evolution (itype — 4, see Section W^ . Such 
tables are thus always dumped on a separate file. 
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call READWT( lun, 'fname', *idniin, *idniax, *nwds, *ierr ) 



Read the weight tables from a disk file. Both the x and /i^ grid must have been defined 
before the call to readwt. On exit the error flag is set as follows: 

Weights are successfully read in. 

1 Read error or input file does not exist. 

2 Input file was written with another QCDNUM version. 

3 Key mismatch (should never occur). 

4 Incompatible x-fi"^ grid definition. 

When successful (ierr = 0), the routine creates the pdf store(s) and returns on exit 
the parameters idmin, idmax and nwds as does the subroutine f illwt. One will get a 
fatal error if there is not enough space in memory to hold all the tables. Like f illwt, 
readwt acts as a do-nothing when a pdf type already resides in memory. 

The code below automatically maintains an up-to-date weight file on disk: 

call readwt (lun, ' polarised. wgt' , idmin, idmax, nw, ierr) 
if (ierr .ne.O) then 

call fillwt(2, idmin, idmax, nw) 

call dmpwgt(2,lun, ' polarised. wgt' ) 
endif 



call NWUSED( *nwtot, *nwuse, *nwtab ) 



Returns the size nwtot of the QCDNUM store (the parameter nwf in qcdnum. inc), the 
number of words used (nwuse) and the size of one pdf table (nwtab). 

5.4 Parameters 

In this section we describe the qcdnum routines to set evolution parameters like the 
perturbative order, flavour thresholds, ttg, etc. All these parameters have reasonable 
defaults but one can change them at any point in the code. Note that a re-definition of 
these evolution parameters invalidates the pdf tables of all existing types. The weight 
tables are not invalidated. In this way, all pdfs are always evolved with the same set 
of parameters; one cannot, for instance, have both un-polarised and polarised pdfs in 
memory, and evolve one in NNLO and the other in NLO. 



call SETORD I GETORD ( lord ) 



Set (or get) the order of the qcdnum calculations to 1, 2 or 3 for LO, NLO and NNLO, 
respectively. Default, lord = 2. 
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call SETALF I GETALF ( alfs, r2 ) 



Set or get for the as evolution the starting value alfs and the starting renormalisation 
scale r2. Default as(^) = 0.118. 



call SETCBT( nfix, iqc, iqb, iqt ) 



nf ix Number of flavours in the FFNS mode. If not set to 3, 4, 5 or 6, QCDNUM 

runs in the VFNS mode. 
iqc,b,t Grid indices of the quark mass thresholds fi^ht- This input is ignored when 

QCDNUM runs in the FFNS mode, that is, when nfix is set to 3, 4, 5 or 6. There 

are some restrictions, dictated by the evolution and interpolation routines: 

iqc > 2, iqb > iqc+2 and iqt > iqb+2. 

A threshold index value of zero (or larger than the number of grid points) means 'beyond 
the upper edge of the grid'. For instance, (iqc,b,t) = (0,0,0) is like running in the 
FFNS with Uf = 3 while the setting (2,4,0) puts the top quark threshold beyond the 
evolution range. By default, QCDNUM runs in the FFNS with Uf = 3. 



call GETCBTC *nfix, *q2c, *q2b, *q2t ) 



Return the current threshold settings. If nfix is non-zero on return, QCDNUM runs in 
the FFNS and the values of q2c,b,t are irrelevant. When nfix = 0, QCDNUM runs in 
the VFNS and the routine returns the threshold values (not the indices) on the /i^ scale. 



call SETABR I GETABR ( ar, br ) 



Define the relation between the factorisation scale /ip and the renormalisation scale fi^ 

2 _ 2 I t 

/^R — CtR/^F + "R- 

Default: ar = 1 and br = 0. 



rscale2 = RFROMF( fscale2 ) fscale2 = FFROMR( rscale2 ) 



Convert the factorisation scale fip to the renormalisation scale fi^^ and vice versa. 
5.5 Evolution 



alphas = ASFUNCC r2, *nf, *ierr ) 



Standalone evolution of ctg on the renormalisation scale yU^ (without using the fi^ grid 
or weight tables). Qcdnum internally keeps track of as so that there is no need to call 
this function; it is just a user interface that gives access to as{fi^). 
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r2 Renormalisation scale /i^ where as is to be calculated. 

nf Returns, on exit, the number of flavours at the scale r2. 

ierr = 1 Too low value of r2. Internally, there is a cut r2 > 0.1 GeV^ and also a 
cut on the slope, to avoid getting too close to A^. 

The input scale and input value of ag, the order of the evolution and the flavour thresh- 
olds are those set by default or by the routines described in Section 15.41 Note that 
although as is evolved on the renormalisation scale the result, in the VFNS, may still 
depend on the relation between fi'^ and /ip. This is because the position of the heavy 
flavour thresholds depends on this relation. 



call EVOLFGC itype, func, def , iqO, *epsi ) 



Evolve a complete set of parton momentum densities from an input scale /Iq. If qcdnum 
runs in the FFNS, the gluon and 2nf quark densities must be given as an input at /Xq. 
In the VFNS, the gluon and 2nf = 6 light quark densities must be given at /Xq < fi"^. 

Here and in the following the parton densities are written on the flavour basis (note the 
PDG convention) with an indexing deflned by 

-6-5-4-3-2-1 1 2 3 4 5 6 

= = (5.1) 

t b c s u dgduscbt 



itype Type of evolution: un-polarised (l), polarised (2), time-like (3), or custom (4). 
func User deflned function func(j,x) (see below) that returns the input parton 

momentum density xfj{x) at iqO. Must be declared external in the calling 

routine. The index j runs from (gluon) to 2nj. 

def Input array dimensioned in the calling routine to def (-6:6,12) which contains 
in def (i, j) the contribution of parton species i to the input distribution j, 
that is, def (i , j ) specifles the flavour decomposition of all input distributions j . 
The indexing of i is given in (15. ip . Internally, qcdnum constructs from def 
a 2nf X 2nf sub-matrix of coefficients and tries to invert this matrix. If that 
fails, the 2nf input densities are not linearly independent in flavour space and 
an error condition is raised, see also (I2.32p . 

iqO Grid index of the starting value /Xq. When evolving in the FFNS the input scale 
can be anywhere inside the range of the /i^ grid. In the VFNS, however, /Iq 
should be below the charm threshold. 

epsi Maximum deviation of the quadratic spline interpolation from linear interpola- 
tion mid-between the grid points (see Section l33|) . A large value epsi > elim 
may indicate spline oscillation and will cause a fatal error message. The value 
of elim can be set by a call to setval. When elim < the error condition 
is disabled so that one can investigate the cause of oscillation. Note that, by 
deflnition, epsi = when QCDNUM is run in the linear interpolation scheme. 

The input function func must be coded as follows 
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double precision function func(ipdf ,x) 
implicit double precision (a-h,o-z) 
if (ipdf .eq.O) then 

func = xgluon(x) 
elseif (ipdf .eq. 1) then 

func = my_f avourite_quark_dstn_l(x) 
elseif (ipdf .eq. 2) then 

func = my_f avourite_quark_dstn_2(x) 
elseif (ipdf .eq. 3) then 



!0 = gluon xg(x) 
! 1 = quarks xql (x) 
!2 = quarks xq2(x) 



endif 

return 

end 

Because evolf g will call func only at the grid points Xi, it is possible to feed tabulated 
values into the evolution routine as is illustrated by the following code 

double precision function pdf input (ipdf ,x) 

implicit double precision (a-h,o-z) 

common /input/ table(0: 12,nxx) Itable with input values 

ix = ixfrmx(x) 

pdf input = table (ipdf ,ix) 

return 

end 



Here is code that evolves both un-polarised and polarised pdfs. 

call fillwtd, idmin, idmax, nw) 
call fillwt(2, idmin, idmax, nw) 

call evolf g(l, fund, defl, iqOl, epsil) 
call evolf g(2, func2, def2, iq02, epsi2) 



lunpolarised 
Ipolarised 

lunpolarised 
Ipolarised 



5.6 External Pdfs 



In QCDNUM, one can read up to 5 different pdf sets from some external source, with type 
identifiers running from 5 to 9. Before reading an external pdf set, care should be taken 
that the perturbative order, the flavour scheme, the positions of the thresholds and the 
input value of as are set correctly in qcdnum. Otherwise one will get the wrong answer 
when the pdf set is used later on in structure function or cross-section calculations. We 
remind that all pdf sets in memory — including the external ones — are invalidated when 
a QCDNUM parameter is re-set by calling one of the routines in Section 15.41 



call PDFINP ( subr, iset, offset, *epsi, *nwds ) 



subr User supplied subroutine (see below), declared external in the calling routine. 
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iset Pdf set identifier in the range 5-9. If tlie pdf set already exists, it will be 

overwritten. 

offset Relative offset at the thresholds fi\. This parameter is used to catch discon- 
tinuities at the thresholds, if any, by sampling the pdfs at /i^(l ± 5). A small 
number like 10~^ should be sufficient, but this depends on how the pdfs are 
externally represented, and how accurate the thresholds are set in qcdnum. 

epsi Maximum deviation of the quadratic spline interpolation from linear interpo- 

lation mid-between the grid points. As for the routine evolsg, a large value 
epsi > elim may indicate spline oscillation and will cause a fatal error mes- 
sage. Note that, by definition, epsi = when qcdnum is run in the linear 
interpolation scheme. 

nwds Last word occupied in the store. Fatal error if the store is not large enough. 

The routine subr provides the interface between qcdnum and the external repository: 

subroutine SUBR ( x, qmu2, xf ) 
implicit double precision (a-h,o-z) 
dimension xf(-6:6) 



The output array xf (-6 : 6) should contain the values of the gluon and the (anti-)quark 
momentum densities at x and /x^, indexed according to (15. ip : note the PDG convention. 



5.7 Pdf Interpolation 

Here we describe routines to access the gluon distribution (xg), the quark and anti-quark 
distributions {xq,xq), or linear combinations of the quarks and anti-quarks. It is also 
possible to directly access the basis singlet /non-singlet pdfs in memory (xe^, defined in 
Section [2.4p . These routines perform local polynomial interpolation on a A; x 3 mesh 
around the interpolation point in x and /x^, where k is the current interpolation order 
in X. Fast routines return the value of a pdf at a given grid point (ix,iq). Two routines 
are provided to investigate the behaviour of the internal spline representation in x. 

In the routines below, the pdf set identifier iset selects the pdf set (or type): un- 
polarised (1), polarised (2), fragmentation function (3), custom (4), or external (5-9). 



Lval = CHKPDFC iset ) 



Returns . true . if the pdf set exists in memory. Both Lval and chkpdf should be 
declared logical in the calling routine. 



pdf = FVALXQ ( iset, id, x, qmu2, ichk ) 



Returns the gluon density or one of the (anti-) quark densities, interpolated to x and /x^. 

iset Pdf set identifier [1-9]. 
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id Gluon, quark or anti-quark identifier, indexed as given in f IS.ip . 

X, qmu2 Input value of x and /i^. 

ichk If set to zero, f valxq will return a null value when x or qmu2 are outside the 

grid boundaries; if set to a non-zero value a fatal error message will be issued. 

The fast version of this function is: pdf = f valij ( iset, id, ix, iq, ichk ). 



call FPDFXQ ( iset, x, qmu2, *pdfs, ichk ) 



Returns all pdf values in one call. The arguments are as given above, except 

pdf s Output array, dimensioned to pdf s (-6 : 6) in the calling routine. The indexing 

is given in (15. ip . 

The fast version is the subroutine f pdf ij ( iset, ix, iq, *pdfs, ichk ). 



pdf = FSUMXQ ( iset, def , x, qmu2, ichk ) 



Return a weighted sum of quark densities. The arguments are as given above, except 

def Input array, dimensioned to def (-6:6) in the calling routine, containing the 

coefficients of the linear combination. The indexing is as given in (15. ip but 
note that def (0) is ignored since it does not correspond to a quark density. 

The fast call is: pdf = f sumij ( iset, def, ix, iq, ichk ). 



pdf = FSNSXQ ( iset, id, x, qmu2, ichk ) 



Return the gluon density or one of the singlet /non-singlet basis pdfs. The arguments 
are as given above, except that id is now indexed as follows: 

0123456789 10 11 12 

; — ; — ; — ; — ; = — : — : — : — - (5-2) 

g Qs 62 63 64 65 Cg q^ 62 63 64 65 Cg 
The fast call is: pdf = f snsij ( iset, id, ix, iq, ichk ). 



pdf = FSPLNE ( iset, id, x, iq ) 



This routine is identical to f snsxq above, except that the local polynomial interpolation 
in X is replaced by spline interpolation, as is done in the qcdnum evolution and convo- 
lution routines (note that f spine does not interpolate in jj?). This function is provided 
as a diagnostic tool to investigate quadratic spline oscillations, if any, which may not 
be visible in the local polynomial interpolation. You do not need this function to detect 
spline oscillations, since that is done automatically by evolsg and pdf inp. 



epsi = SPLCHK ( iset, id, iq ) 



Returns e = ||w — t;|| at a grid point iq. Here u and v are the vectors of quadratic and 
linear interpolation mid-between the grid points in x, as is described in Section [331 By 
definition, e = for linear interpolation, and should be a small number (like 0.05, say) 
for quadratic interpolation. Large values indicate that the spline oscillates. 
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6 Convolution Engine 

The QCDNUM convolution engine provides tools to calculate structure functions in deep 
inelastic scattering, hadron-hadron scattering cross-sections and parton luminosities. 
The engine drives the add-on package zmstf that computes the zero-mass structure 
functions F2, F]^ and xF^, in un-polarised deep inelastic scattering. It is also used in 
the HQSTF package that computes the heavy flavour contributions to F2 and Fl in the 
fixed flavour number scheme [1^. Both these packages are included in the qcdnum 
distribution and are described in the Sections IC.31 and ID. II of this write-up. 

From the parton number densities / and kernels K, all kind of convolution integrals can 
be calculated with the engine, such as 

x[f®K]{x), x[f®Ka®Kj,]{x), X[fa®h]{x), x[f a ® fb ® K]{x) , CtC. 

Here ® stands for Mellin convolution as defined by (12.51) . We refer to Section [3l2l for 
how convolution integrals are computed and where the factor x in front comes from. We 
emphasise that the kernel K must be defined by convolution with a number density. If 
not, then it must be transformed as necessary, before it is fed into qcdnum. 

The steps to be taken in a calculation based on the convolution engine are the following. 

1. Declare one or more stores and partition these into tables. Then fill the tables 
with weights for all the convolution kernels needed in the calculation (Section 16.21) : 



2. Write a function myfun(ix,iq) that returns the structure function, cross section 
or luminosity at a grid point in x and /i^ (Section l6.3p : 

3. Pass myf un to a QCDNUM routine that will take care of the interpolation to any 
desired x and /x^ (Section 16.41) . 

This procedure is fairly straight-forward and therefore suitable for prototyping and de- 
bugging. However, there is a considerable amount of overhead so that it is recommended 
to ultimately move steps 2 and 3 of the computation to a fast calculation scheme that is 
described in Section 16. 5[ By this one will gain at least an order of magnitude in speed. 

Before we present the convolution engine we will first, in the next section, introduce the 
rescaling variable x to accommodate generalised mass variable flavour number schemes 
(gm-vfns, see [27] for a recent review) in structure function calculations. 

6.1 Rescaling Variable in Convolution Integrals 



The general expression for a structure function can be written as 

'X 



J^i{x,Q'^) = ^x —fj{z,fx^)C^j -,/x^Q^m^,as(/^^) • (6.1) 

■ Jy Z LZ J 



Here the index i labels the structure function {e.g. F2, Fl, xF^, Fg, . . .) and j labels a 
parton number density like the gluon, the singlet and various non-singlets. The coeffi- 
cient function Cij depends on x, on the scale variables //^ and Q^, on one or more quark 
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masses m^ and on the strong coupling constant as- The variable x = ^2;, a > 1, is a 
so-called resettling variable which takes into account the kinematic constraints of heavy 
quark production, for instance, 

X = ax=(l + ^^ X. (6.2) 

We have < x < 1 so that the range of x in (16. ip is restricted to < x < 1/a. In the 
zero- mass limit a = 1, x = 2;, and (16. ip reduces to the Mellin form x[/ ® C]{x). 

To calculate the structure function, we first have to evaluate the convolution integrals 
(for clarity we drop a^ and the indices i,j) 



As in Section 13.21 we denote by h{y, t) a parton momentum density in the logarithmic 
scaling variables y = —In a; and t = In/x^. In terms of these, and provided that x is 
proportional to x, (16. 3p can be written as a weighted sum of spline coefficients 

i 

J'{y^,Q') = Y,W^J^^ (6-4) 

with Wij = Wi-j+i and 

ryt—b 
we = e-^ dzYi{z)D{ye-b-z,t,Q^,ml) (1 < £ < n). (6.5) 

Jo 

Here D{y,t,Q'^,m\) = e~^C(e~^, e*, (5^,m^) and h = ln(a). It is understood that the 
integral (16. 5 p is set to zero in case ye — b < 0. In the massive schemes, b > depends 
on t which implies that the weights must be stored in 2-dimensional y-t tables. 

We emphasise that convolution integrals found in the literature must, if necessary, be 
brought into the general form (16. ip by modifying the published Wilson coefficient. An 
example of such a modification can be found in Appendix [Dl 



6.2 Weight Tables 

In this section we describe routines that partition a linear store into tables and fill these 
tables with weights used in the calculation of convolution integrals. It is important to 
realise that the convolution kernels may contain singularities, see also Appendix |Al To 
deal with such singularities, we formally decompose a kernel into a regular part (A), a 
singular part (B), a product (RS) and a delta function 

C{x) = A{x) + [B{x)]+ + R{x)[S{x)]+ + D{x)6{l - x). (6.6) 

QcDNUM provides routines that can calculate weights for each term separately (if 
present) and add these to the weight table of C. 

For reasons of efficiency and economy of storage, there are four different types of tables: 
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itype 
itype 
itype 
itype 



1 
2 

3 

4 



Weights that depend only on x. Table identifiers run from 101-199; 
Weights that depend on x and Uf. Identifiers run from 201-299; 
Weights that depend on x and /x^. Identifiers run from 301-399; 
Weights that depend on x, /i^ and nf. Identifiers run from 401-499. 



Although it is a good idea to take out as many /x^-dependent factors as possible from the 
convolution kernel, it is clear from fl6.3l) that quark mass parameters and the relation 
between /i^ and Q^ may enter via the rescaling variable x ^i^nd that this dependence can 
never be factored out of the convolution integral. Thus the weight tables of the GM 
schemes will, in general, depend on x and /x^ and must be stored in type-3 or 4 tables. 

QCDNUM calculates by Gauss quadrature (cernlib routine D103) the integrals that 
define the weights. In case the default accuracy of e = 10~^ cannot be reached (fatal error 
message), this limit can be raised by a call to setvaK'epsg' , value). Note, however, 
that problems with the Gauss integration will most likely be caused by problems with the 
integrand — such as near-singular behaviour somewhere in the integration domain — and 
that this cannot be cured by relaxing the required accuracy. 

In Table [3] we list all available weight routines. 

Table 3: Qcdnum convolution weight table routines. 



Subroutine or 


function 


Description 


BOOKTAB ( 


; w, 


nw, itypes, *nwords ) 


Partition into tables 


MAKEWTA ( 


'- w, 


id, afun, achi ) 


Regular piece A{x) 


MAKEWTB ( 


; w, 


id, bfun, achi, nodelta ) 


Singular piece [B{x)] + 


MAKEWRS ( 


; w, 


id, rfun, sfun, achi, nodelta ) 


Product R{x)[S{x)] + 


MAKEWTD ( 


; w, 


id, dfun, achi ) 


Delta function D{x)6{l — x) 


MAKEWTX ( 


; w, 


id ) 


Weight table for x[fa ® fb 


SCALEWT ( 


'- w. 


c, id ) 


Scale weight table 


IDSPFUN ( 


: 'p 


ij ' , lord, itype ) 


Splitting function index 


COPYWGT ( 


; w, 


idl, id2, iadd ) 


Copy weight table 


WCROSSW ( 


; w, 


ida, idb, idc, iadd ) 


Double convolution weights 


WTIMESF ( 


; w, 


fun, idl, id2, iadd ) 


Multiply by /(/i2,n^) 


SETWPAR ( 


; w. 


pars , n ) 


Store extra information 


GETWPAR ( 


'- w. 


*pars, n ) 


Read extra information 


TABDUMP ( 


'- w, 


lun, 'filename', 'key' ) 


Dump to disk 


TABREAD ( 


'- w, 


n, lun, 'fn', 'key', *nw, *ierr 


) Read from disk 



Output arguments are pre-fixed with an asterisk (*). 



call BOOKTAB ( w, nw, itypes, *nwords ) 



Partition a store w into tables. 
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w Double precision array declared in the calling routine. 

nw Dimension of w as declared in the calling routine. 

itypes Integer array dimensioned to itypes(4) in the calling routine which contains 

in itypes (i) the number of tables (< 99) of type i to be generated. When 

itypes (i) = then no tables of type i will be generated. 
nwords Gives, on exit, the number of words used in the store. If nwords is negative, 

then the store is not sufficiently large and should be re-dimensioned in the 

calling routine to at least -nwords. 

Note that one can declare and partition as many stores as desired, one per structure 
function for instance. 



call MAKEWTA ( w, id, afun, achi ) 



Calculate the weights for the regular contribution A{x) to a convolution kernel and add 
these to table id in the store w. 



w Store declared in the calling routine and previously partitioned by booktab. 

id Table identifier. To add results to a type-n table, one should use identifiers in 

the range n01-n99, with n = 1, 2, 3 or 4. 
afun User function (see below) returning the regular piece of the convolution kernel. 

Should be declared external in the calling routine, 
achi User function (see below), declared external in the calling routine, that returns 

the value a of the rescaling variable x = (^^■ 

The function afun provides an interface between qcdnum and the regular part of the 
kernel C(x,/i^, Q^,""^^) and should be coded as follows|if| 

double precision function afun(chi,qmu2,nf ) ! chi = a*x 
implicit double precision (a-h,o-z) 

common /fixpar/ pari, par2, 

Q2 = some_function_of (qmu2,some_params) 
afun = cfun(chi,qmu2,Q2,nf ,some_params) 
return 
end 



! parameters, if any 

!Q2 

! convolution kernel 



The function achi should return, as a function of /i^, the factor a that defines the 
rescaling variable x = (^^■ 

double precision function achi(qmu2) 
implicit double precision (a-h,o-z) 

common /fixpar/ pari, par2, [parameters, if any 

Q2 = some_function_of (qmu2,some_params) ! Q2 
achi = some_function_of (Q2,some_params) 
return 
end 
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We assume here that the kernel conforms to (j6.ip . If not, then afun must take care of this. 



45 



QcDNUM insists that always achi > 1, one will get a fatal error if not. To compute 
standard Mellin convolutions x[/ (S) C]{x), simply set achi = 1 for all /i^. 

double precision function achi(qmu2) 

implicit double precision (a-h,o-z) 

achi = l.DO 

return 

end 



call MAKEWTB ( w, id, bfun, achi, nodelta ) 



Calculate the weights for the singular contribution [i?(x)]+ to a convolution kernel and 
add these to a table in the store w. The arguments and the coding of bfun and achi 
are as for makewta. Thus, if a kernel has both a regular and a singular part, then do 

call makewa(w, 201 ,afun, achi) !put weights in id = 201 
call makewb ( w, 20 1 , bfun, achi, 0) !add weights to id = 201 

It is seen from Appendix |Al equation (JA.41) . that a '+' prescription generates a 6{1 — x) 
contribution. By default, makewtb includes this contribution, unless you set nodelta = 1. 
In that case the S{1 — x) contribution is not calculated and must be entered, perhaps 
combined with other such contributions, via a call to makewtd, see below. 



call MAKEWRS ( w, id, rfun, sfun, achi, nodelta ) 



Calculate the weights for the product contribution i?(x)[S'(x)]+ to a convolution kernel 
and add these to a table in the store w. The arguments and the coding of rfun, sfun 
and achi are as for makewta. 



call MAKEWTD ( w, id, dfun, achi ) 



Calculate the weights for the S{l—x) contribution to a convolution kernel and add these 
to a table in the store w. The delta function is multiplied by the function dfun. The 
arguments and the coding of dfun and achi is as for makewta. 



call MAKEWTX ( w, id ) 



Calculate the weights f l3.20p for the convolution x[fa ® fb]{x). 

w Store declared in the calling routine and previously partitioned by booktab. 

id Table identifier. Because the weight table depends only on x, it can be stored 

in a type-1 table, but equally well in types-2, 3 or 4, if desired. 
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call SCALEWT ( w, c, id ) 



Multiply the contents of table id by a constant c. 



id = IDSPFUN ( 'pij', lord, itype ) 



Return the index (< 0) of a splitting function weight table stored internally in qcdnum. 
'pij ' Name of the splitting function. Valid input strings are 
PQQ, PQG, PGQ, PGG, PPL, PMl, PVA. 



lord Select LO (1), NLO (2) or NNLO (3). 

itype Select evolution type: un-polarised (1), polarised (2), fragmentation func- 
tion (3) or custom (4). 

The index returned by idspfun is encoded as -(1000*itype+id), where id is the 
internal table identifier. If the table does not exist, the function returns a value of -1. 



call COPYWGT ( w, idl, id2, iadd ) 



Copy the contents of table idl to id2. 

w Store declared in the calling routine. 

idl Input table identifier. One can copy a splitting function weight table from 

internal QCDNUM memory to the store by setting idl < 0. Valid identifiers are 

generated by idspfun (), as described above. 

id2 Output table identifier with id2 7^ idl. The output table type may be different 

from the input table type, see below, 
iadd If set to copy idl to id2, if set to +1 (-1) add (subtract) idl to (from) id2. 

For this routine — and for those described below — the output table type can be different 
from the input table type, provided that this does not lead to a loss of input information. 
Thus one can copy a type-1 table to a type-3 table but not the other way around (fatal 
error). Note that input splitting function weight tables are all type-2. 



call WCROSSW ( w, ida, idb, idc, iadd ) 



This routine generates a weight table for the convolution of two kernels Ka and Ki,. The 
weight table is calculated with 03.18P from two input tables Wa and W^,. 

w Store declared in the calling routine. 
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ida Table identifier containing the weights of kernel Ka- When ida < one will 
access a splitting function weight table which is stored internally in qcdnum. 
See idspfunO above for how to generate a valid splitting function identifier. 

idb As above for the weights of kernel Ki,. 

idc Output table identifier. Cannot be set equal to ida or idb. 

iadd If set to store the result of the convolution in idc, if set to +1 (-1) add 
(subtract) the result to (from) the contents of idc. 

The table types of ida and idb may be different, but the type of idc must be such that 
it can contain either input table. Thus if ida is type-2 (x, n/) and idb is type-3 (x, /x^), 
then idc must be type-4 (x, /i^,n/). The routine checks this. 



call WTIMESF ( w, fun, idl, id2, iadd ) 



Multiply a weight table by a function of /i^ and nf and store the result in another table. 

w Store declared in the calling routine. 

fun User supplied double precision function fun(iq,nf ) declared external in the 
calling routine. 

idl Input weight table identifier. It is possible to access a splitting function weight 
table by setting idl < 0. Valid identifiers can be obtained from idspfunO 
described above. 

id2 Identifier of the output table. It is allowed to have idl = id2 (in-place modi- 
fication of a table), unless idl is a splitting function table. The table type of 
id2 must be such that no information is lost. The routine checks this. 

iadd Store the result in id2 in case iadd = or add (subtract) the result to (from) 
id2 in case iadd = +1 (-1). 

The routine loops over iq and nf and calls fun(iq,nf ) with the following argument 
ranges, depending on the output table type: 



type 


variables 


iq range 


nf range 


1 


X 




1-1 


3-3 


2 


X, Uf 




1-1 


3-6 


3 


X, /i^ 




1-nq 


3-3 


4 


X, /i^. 


^/ 


1-nq 


3-6 



With this routine one can, in combination with wcrossw, construct weight tables for 
combinations of convolution kernels, such as those given in (lC.7p . For instance, here is 
code that generates a table for 

external betaO !beta function 

call WcrossW ( w, idC2Q0, idSpfun('PPL' ,2, 1) , idC2P21, ) 
call WcrossW ( w, idC2Pl, idSpfun('PQQ' , 1, 1) , idC2P21, +1 ) 
call WtimesF ( w, betaO , idC2Pl , idC2P21, -1 ) 
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call SETWPAR ( w, par, n ) 



Write extra information to the store, for instance quark masses or other parameters that 
you may want to dump to disk, together with the tables themselves. 

w Store, partitioned by a previous call to booktab. 

par List of parameters to be written. Should be dimensioned to at least par (n) 

in the calling routine, 
n Number of items to be written up to a maximum of miwO = 20. If necessary, 

one can change the value of miwO in qcdnum. inc and recompile QCDNUM. 

The parameters can be read back by a call to getwpar(w,par,n). 



call TABDUMP ( w, lun, 'filename', 'key' ) 



Dump the store w to disk. Apart from the store, information is written about the 
QCDNUM version, the x-/i^ grid definition and the current spline interpolation order. The 
key text string can be used to stamp the file with a version number or other identifier. 
The dump is unformatted so that the file cannot be exchanged across machines. 



call TABREAD ( w, nw, lun, 'filename', 'key', *nwords, *ierr ) 



Read a store from disk into the array w(nw). The size of the store (in words) is returned 
in nwords. You will get a fatal error message if w(nw) is not large enough to contain the 
store. Note that the x and /i^ grids must have been defined before the call to tabread. 
On exit, the error fiag is set as follows (non-zero means that nothing has been read in). 

Store successfully read in. 

1 Read error or input file does not exist. 

2 File written by another qcdnum version. 

3 Key mismatch. 

4 Incompatible x-^"^ grid definition. 

Qcdnum insists that the key written on the file matches the key entered as an argument 
to tabreadO Thus if, for instance, the key is set to a package name and version number 
then the user of the package cannot read obsolete files written by earlier versions, or 
read files written by another package. If you don't want to use keys, just enter an empty 
string as a key in the calls to tabdump and tabread. 



-'^^Note that the key matching is case insensitive and that leading and trailing blanks are ignored. 
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6.3 Convolution 

In Table H] we list the routines that can be used to build a structure function, cross- 
Table 4: Calls in the QCDNUM convolution engine. 

Subroutine or function Description 

FCROSSK ( w, idw, iset, idf, ix, iq ) Convolution x[/ ® -ft'] 

FCROSSF ( w, idw, iset, ida, idb, ix, iq ) Convolution x[/a ® /h] 

EFROMQQ ( qvec, *evec) Transform from g, g to e^ 

QQFROME ( evec, *qvec) Transform from e^ to g, g 

NFLAVOR ( iq ) Returns Uf 

GETALFN ( iq, n, *ierr) Returns (as/27r)" 

Output arguments are pre-fixed with an asterisk (*). 

section or parton luminosity at a grid point in x and /i^. 

A convolution is always computed with the pdfs in qcdnum memory, that is, with the 
gluon density or with one of the singlet/non-singlet quark densities |e^) as defined in 
Section 12.41 To translate a linear combination of quarks and anti-quarks to the |e^) 
basis, and vice versa, the routines ef romqq and qqf rome are provided. 

For convenience we show here again the indexing (15.21) of the singlet/non-singlet basis 



0123456789 10 11 12 
g Qs 62 63 64 65 Cg q^ 62 63 64 65 eg 

and the indexing (15.11) of the flavour basis 

-6-5-4-3-2-1 1 2 3 4 5 6 
i b c s u dgduscht 



(6.7) 



(6.8) 



val = FCROSSK ( w, idw, iset, idf, ix, iq ) 



Calculate the convolution x[f ® K]{x) at a grid point in x and ^j?. 

w Store declared in the calling routine and previously filled with weights. 

idw Identifier of a table in the store w. 

iset Pdf set identifier [1-9] 

idf Pdf identifier, indexed according to (16.71) . 

ix, iq Indices of an x-//^ grid point. 

Splitting function tables cannot be directly accessed by this routine; they should first 
be copied to the store by a call to copywgt. 



^^ Un-polarised (1), polarised (2), fragmentation function (3), custom (4), or external (5-9). 

50 



val = FCROSSF ( w, idw, iset, ida, idb, ix, iq ) 



Calculate the convolution x[fa CS) fb]{x) at a grid point in x and /x^. 

w Store declared in the calling routine and previously partitioned by booktab. 

idw Identifier of a weight table, previously filled by a call to makewx. 

iset Pdf set identifier [1-9]. 

ida, idb Pdf identifiers, indexed according to (16.71) . 

ix, iq Indices of an a;-/i^ grid point. 

A convolution of a linear combination of pdfs must be calculated as a sum of pair-wise 
convolutions, with each term computed by f crossf . Note that this is much easier done 
with the fast routines described in Section 16.51 



call EFROMQQ ( qvec, *evec, nf ) 



Transform the coefficients of a linear combination of quarks and anti-quarks from the 
flavour basis to the singlet /non-singlet basis as described in Section 12.41 

qvec Input array, dimensioned qvec (-6: 6), filled with the coefficients of a linear 
combination of quarks and anti-quarks and indexed according to fl6.8p . 

evec Output array, dimensioned to evec(12), filled with the coefficients written on 
the singlet /non-singlet basis, indexed according to (16. 7p . 

nf Active number of fiavours. This parameter is needed to construct the appropri- 

ate 2n/ X 2nf transformation matrix that acts on qvec(-nf :nf ). 

Thus if a linear combination of quarks and anti-quarks is written as 

\p) = J2(»^h) + m)) = J2('^t\et) + dr\er)) (6.9) 

2=1 i=l 

and tti and /3j are stored in the input vector qvec, then the coefficients df are returned 
in the output vector evec. 



call QQFROME ( evec, *qvec, nf ) 



Transform the coefficients of a linear combination of basis vectors from the singlet/non- 
singlet basis to the fiavour basis. The arguments are as for ef romqq. 



nf = NFLAVOR ( iq ) 



Returns the number of active fiavours at the grid point iq. Note that this number 
is (4,5,6) and not (3,4,5) at the thresholds (iqc,iqb,iqt). 
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as = GETALFN ( iq, n, *ierr ) 



Returns the value of (as/27r)" at the factorisation scale ftp. Here ctg is computed from the 
Taylor expansion (12.171) . appropriately truncated depending on the current value of the 
perturbative order (lord). Because the truncation is different for the pdfs (Section 12.31) 
and the structure functions (Appendix IC.2p . the value of n must be set as follows. 



• If the convolution at (LO, NLO, NNLO) should be multiplied by (as, a'^, a^) then 
you set n = (1,2,3) in the call to getalfn. 

• If the convolution at (LO, NLO, NNLO) should be multiplied by (l,as5«s) then 
you set n = (0,-1,-2) in the call to getalfn. 

• If n > lord, the value of (as/27r)" is calculated at /z^, instead of at /x|. (This is 
already the case for n = lord, see Section 1^751 ) 

To have access to the NNLO discontinuities at the thresholds, one can set iq positive 
(includes discontinuity) or negative (does not include discontinuity), thus: 

call getalfn ( iqcharm, n, ierr ) ! result for nf = 4 

call getalfn ( -iqcharm, n, ierr) ! result for nf = 3 

In other words, by preceding iq with a minus sign one effectively changes the qcdnum 
default Uf = (4, 5, 6) at the thresholds to the alternative nj = (3, 4, 5). When iq is close 
to or below the value of A^, then ierr = 1 and getalfn returns the null value. This 
also happens when iq is outside the grid boundaries (ierr = 2). 

Note that the qcdnum expansion parameter is as/27r but that many convolution kernels 
found in the literature are defined for an expansion in ag/^n, in which case one must 
account for the appropriate factors of 2" somewhere in the calculation. 



6.4 Interpolation 

With the routines presented above one can write a function stfun(ix,iq) that returns 
the value of a structure function or cross section at a grid point in x and /i^. The routine 
stf unxq then takes care of the interpolation to any value of x and /x^. This interpolation 
is done on a fc x 3 mesh around the interpolation point, where k is set to the current 
spline interpolation order (2 = linear, 3 = quadratic). Thus 3k functions have to be 
computed for each interpolation which becomes inefficient if there are interpolations 
with overlapping meshes. By processing lists of interpolation points, instead of each 
point individually, redundant calculations are avoided which can lead to considerable 
gains in computing time. In other words, stf unxq should not be called in a loop over 
interpolation points but be given the list of points. 



call STFUNXQ ( stfun, x, qmu2, stf, n, ichk ) 



Interpolate the function stfun(ix,iq) to a list of x and /i^ values. 
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stfun Double precision function, declared external in the calling routine, that 

returns the value of a structure function or cross-section at (ix,iq). 

X, qmu2 List of interpolation points, dimensioned to at least n in the calling routine. 

stf Contains, on exit, the list of interpolated results. 

n Number of items in x, qmu2 and stf. 

ichk If set to the routine returns a null value if x or /x^ are outside the 

boundaries of the grid; if set non-zero it will insist that all interpolation 
points are inside the grid boundaries. 

Note that the interpolation is done in //^ and not in Q^. You have to keep track yourself 
of the relation between these two scales. 



6.5 Fast Computation 

As already remarked above, the routines provided up to now are fine for prototyping but 
are slow because there is quite a lot of overhead when the calculation is repeated at more 
than one interpolation point. Here we describe a set of routines that does optimised 
bulk calculations on selected points in the x-fi"^ grid. With these fast routines one can 
easily gain one or two orders of magnitude in speed. To make the calculation flexible 
it is broken down into small steps where intermediate results are stored into scratch 
buffers (by default, the fast engine generates 5 scratch buffers but one can have more, 
if necessary). An optimised calculation proceeds as follows. 

1. Pass a list of interpolation points in x and yU^ to a qcdnum routine that determines 
which grid points will be occupied in the course of the calculation; 

2. Store a pdf or a linear combination of pdfs in a scratch buffer; 

3. Convolve the pdf with a convolution kernel or a perturbative expansion of kernels; 

4. Multiply the convolution by a function of x and /x^, for instance by some power 
of Os or by some kinematic factor; 

5. Accumulate the cross section or structure function in a final buffer; 

6. Pass this buffer to an interpolation routine to get a list of interpolated results. 

The list of subroutines is given in Table [51 In principle, the output buffer of any fast 
routine can serve as the input buffer of any other fast routine. There is, however, a 
little complication related to the amount of information stored in a buffer. For interpo- 
lation purposes, it is sufficient to store results only at the mesh points; such a buffer is 
called sparse. A convolution routine, on the other hand, does not only need the values at 
the mesh points Xj, but also the values at all points Xj > Xj. An input buffer with such 
a storage pattern is called dense; a dense buffer is of course more expensive to generate 
than a sparse buffer. Usually one does not have to worry about sparse and dense buffers, 
because qcdnum has reasonable defaults on what kind of buffer is accepted as input, 
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Table 5: Fast convolution engine. 



Subroutine or function 


Description 


FASTINI ( 


[ X, qmu2, n, ichk ) 


Pass list of X and /i^ values 


FASTCLR ( 


: id ) 


Clear buffer 


FASTEPM ( 


[ iset, idf, idout ) 


Store \g, e^) in a scratch buffer 


FASTSNS ( 


[ iset, pdf, isel, idout ) 


Store singlet /non-singlet component 


FASTSUM ( 


; iset, coef, idout ) 


Store weighted sum of e^) 


FASTFXK ( 


; w, idw, idf, idout ) 


Convolution x[f ^ K]{x) 


FASTFXF ( 


; w, idw, ida, idb, idout ) 


Convolution x[fa ^ fb]ix) 


FASTKIN ( 


; id, fun ) 


Scale by a kinematic factor 


FASTCPY ( 


; idin, idout, iadd ) 


Copy or accumulate result 


FASTFXQ ( 


; id, *f, n ) 


Interpolation 



Output arguments are pre-fixed with an asterisk (*). 



and what kind of buffer is generated on output. One can aways override the output 
default and force a routine to generate a dense or a sparse buffer, as needed. 

The pdf set parameter iset in the routines f astepm, f astsns and f astsum selects the 
pdf set, namely, un-polarised (1), polarised (2), fragmentation function (3), custom (4), 
or external (5-9). 



call FASTINI ( x, qniu2, n, ichk ) 



Pass a list of interpolation points and, at the first call, generate the set of scratch buffers. 



X 

qmu2 

n 

ichk 



Array, dimensioned to at least n in the calling routine, filled with x values. 

As above, but for fi^ (not Q^). 

Number of entries in x and qmu2. 

If non-zero, f astini insists that all x and fi^ are within the grid boundaries. 



By default, 5 scratch buffers (id = 1-5) are generated at the first call (or cleared if they 
exist). This number can be changed by calling setintC'ntab' ,ival) prior to fastini. 
One will get a fatal error if there is not enough space for the scratch buffers, in which 
case one has to increase the value of nwfO in qcdnum.inc, and recompile QCDNUM. 



call FASTCLR ( id ) 



Clear a scratch buffer. Setting id = will clear all buffers. 



call FASTEPM ( iset, idf, idout ) 



Copy the gluon density or one of the basis pdfs |e^) to a scratch table. 
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iset Input pdf set identifier [1-9]. 

idf Pdf identifier [0-12], indexed according to (16. 7p . 

idout Output scratch table identifier [1-5]. 

By default, f astepm generates a dense buffer; a sparse buffer is generated when you 
pre-pend the output identifier with a minus sign. 

call fastEpmd, 0, 1) Idense table output 

call fastEpmd, 0, -1) ! sparse table output 



call FASTSNS ( iset, pdf, isel, idout ) 



Decompose a given linear combination of quarks and anti-quarks into singlet and non- 
singlet components and copy a specific component to a scratch buffer. 

iset Input pdf set identifier [1-9]. 

pdf Input array, dimensioned pdf (-6:6), filled with the coefficients of a linear 

combination of quarks and anti-quarks and indexed according to fl6.8p . 

isel Selection flag [0-7], see below. 

idout Output scratch table identifier [1-5]. 

The isel flag selects the gluon density (0), the singlet component gs (1), the non-singlet 
component g+ (2), the valence component q^ (3), the non-singlet component g~g (4), 
the sum q^ + q^^ (5), all non-singlets q^ + q'^^ + g^^ (6) or all quarks (7). 

Note that the singlet is weighted by an appropriate rtf dependent factor, for instance by 
the average square of the quark charges in case pdf corresponds to the charge weighted 
sum of quarks and anti-quarks. Note also that the gluon density is multiplied by the 
same factor, as is required in structure function calculations, see Appendix O Setting 
isel = or 1 is thus not the same as calling f astepm for the identifiers or 1. 

By default, f astsns generates a dense buffer; a sparse buffer is generated when you 
pre-pend the output identifier with a minus sign. 



call FASTSUM ( iset, coef , idout ) 



Copy a linear combination of basis pdfs le"*^) to a scratch table. 

iset Input pdf set identifier [1-9]. 

coef Array of coefficients dimensioned coef (0: 12,3:6) in the calling routine. 

idout Output scratch table identifier [1-5]. 

The array coef (i,nf) is indexed according to fl6.7p . Here is code that fills coef by 
transforming a set of quark coefficients from flavour space to singlet /non-singlet space: 
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dimension qvec(-6:6), coef(0: 12,3:6) 

do nf = 3,6 

coef(0,nf) = O.DO Igluon coefficient 

call efromqqCqvec, coef(l,nf), nf) ! quark coefficients 

enddo 

By masking out coefficients, one can copy the singlet component, or various combinations 
of non-singlets; this is exactly what f astsns does. To copy the gluon distribution, one 
must set all coefficients to zero, except coef (0,nf ). 

By default, fastsum generates a dense buffer; a sparse buffer is generated when you 
pre-pend the output identifier with a minus sign. 



call FASTFXK ( w, idw, idf, idout ) 



Calculate the convolution x[f ® K]{x) at all selected grid points. 

w Store, declared in the calling routine and previously ffiled with weights. 

idw Set of weight identifiers, declared idw (4) in the calling routine, see below, 

idf Input scratch buffer, previously ffiled by f astepm, f astsns or fastsum. 

idout Output scratch table with idout 7^ idf. 

One can either convolve with a given weight table or with a perturbative expansion of 
weight tables, depending on what one puts in the array idw: 

1. To convolve with a given weight table, set idw(l) to the identifier of that weight 
table and set idw (2), idw (3) and idw (4) to zero; 

2. To convolve with a perturbative expansion, store the (LO,NLO,NNLO) weight 
table identifiers in idw(l), idw(2) and idw(3). Set the identifier to zero if no 
such table exists, as is the case for Fl at LO, for instance. Declare in idw (4) the 
leading power of Og, that is, multiply (LO,NLO,NNLO) by (1, as, a^) if idw (4) = 
and by (og, «s5 c^s) i^ idw (4) = 1. Note that the perturbative expansion is summed 
up to the current perturbative order, as defined by an upstream call to setord. 

The routine only accepts a dense buffer as input (otherwise fatal error) and will, by 
default, generate a sparse buffer as output. If you pre-pend the output identifier with a 
minus sign, the output buffer will be dense. In this way, the output table can serve as 
an input to another convolution which allows one to calculate multiple convolutions in 
a chain. For example o 



call fastSumC 1, coef, 1 ) 
call fastFxKC w, idKl, 1, -2 ) 
call fastPxKC w, idK2, 2, 3 ) 



1 = f 

2 = f * Kl 

3 = f * Kl * K2 



^^It is more efficient, however, to first calculate with wcrossw a weight table for K^ — Ki (g) K2, and 
use that table to convolve K3 with /. 
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call FASTFXF ( w, idx, Ida, idb, idout ) 



Calculate the convolution x[fa ® fb]{x) at all selected grid points. 

w Store, declared in the calling routine and previously partitioned by booktab. 

idx Identifier of a weight table, previously filled by a call to makewtx. 

ida, idb Identifiers of input scratch tables. It is allowed to have ida = idb. 

idout Output scratch table with idout ^ ida or idb. 

As above, the routine accepts only dense buffers as input, and generates a sparse buffer 
as output, unless the output identifier is pre-pended by a minus sign, thus. 



call fastSumC 1, coefa, 1 ) 

call fastSumC 1, coefb, 2 ) 

call fastPxFC w, idwX, 1, 2, -3 ) 

call fastFxKC w, idwK, 3, 4 ) 



1 = fa 

2 = fb 

3 = fa * fb 

4 = fa * fb * K 



call FASTKIN ( id, fun ) 



Multiply the contents of a scratch table by a kinematic factor. 

id Identifier of the input scratch table. 

fun Double precision function, declared external in the calling routine. 

The routine loops over the selected grid points and calls the user supplied function fun 
that should return the kinematic factor. The syntax of fun is 

double precision function fun ( ix, iq, nf, ithresh ) 

ix,iq Grid point indices. 

nf Number of fiavors at iq. This number is bi-valued at the thresholds so that 

at the charm threshold, for instance, nf can be either 3 or 4o 

ithresh Set to if iq is not at a threshold and to +1 (-1) if iq is at a threshold 
with the upper (lower) number of fiavours. This variable can be used to 
take NNLO discontinuities into account, as is shown in the example below. 

double precision function fkin(ix,iq,nf , ithresh) 

if (ithresh. ge.O) then 

alfas = getalfnC iq,l,ierr) !alfas/2pi with discontinuity 
else 

alfas = getalfn(-iq, Ijierr) ! without discontinuity 
endif 



22 



The reader may wonder when QCDNUM returns the value 3, and when the value 4. This depends 



on the interpolation point /i^ to which iq is associated: if //^ is below (above) /i^, then nf = 3 (4). 
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call FASTCPY ( idin, idout, iadd ) 



Copy or accumulate a result in an output buffer. 

idin Identifier of the input scratch table. 

idout Identifier of the output scratch table with idout ^ idin. 

iadd Store (O), add (l) or subtract (-1) the result to idout. 

The type of output buffer (sparse or dense) is the same as that of the input buffer, 
except that once you have used a sparse input buffer, the output buffer will be flagged a 
sparse and will remain so until you set iadd = to start a new accumulation in idout. 



call FASTFXQ ( id, *f, n ) 



Interpolate the contents of id to the list of x and /i^ values that was previously passed 
to QCDNUM by the call to f astini. 

id Identifier of an input scratch buffer. 

f Array dimensioned to at least n in the calling routine that will contain, on 

exit, the interpolated values. 

n Number of interpolations requested. 

The routine works through the list of interpolation points given in the call to f astini 
and exits when it reaches the end of that list or when the number of interpolations is 
equal to n, whatever happens first. A dense input buffer is allowed, but wasteful since 
it contains a lot of information that is not used by f astf xq. 

6.6 Custom Evolution 

In Section 15.31 we have described the f illwt routine to generate weight tables (and 
pdf tables) for the evolution of un-polarised pdfs (itype = 1), polarised pdfs (2), and 
fragmentation functions (3). But qcdnum can handle yet another type of evolution, 
with user-defined evolution kernels (custom evolution, itype = 4). For this one has to 
provide a subroutine, described below, that generates the weight tables. This routine 
is passed to qcdnum as an argument of the custom weight filling routine f illwc, after 
which the custom evolution becomes available by switching to itype = 4. 

external mysub, func 

call f illwc ( mysub, idmin, idmax, nwords ) 

call evolfgC 4, func, def , iqO, epsi ) 
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subroutine myweight (w,nw,nwords, idpij ,mxord,idum) 



C-1 
C-2 



C-3 

C-4 



C-5 



w (in) qcdnum store passed by reference 

nw (in) number of words available 

nwords (out) number of words used < not enough space 

idpij (out) list of Pij table identifiers 

mxord (out) maximum perturbative order LO,NLO,NNLO 

idum not used at present 

implicit double precision (a-h,o-z) 

dimension w(*) , idpij (7,3), itypes(4) 

external AChi 

external PQQR, PQQS, FQQD ! PQQ 

external PQGA, PGQA ! PQG, PGQ 

external PGGA, PGGR, PGGS, PGGD ! PGG 

!s/r name for error messages 



call setUmsg( 'myweight ' ) 

Max perturbative order 

mxord = 1 

Partition 

itypes(l) = 2 

itypes(2) = 2 

itypesO) = 

itypes(4) = 

call BookTab(w,nw, itypes , nwords) 

Not enough space 

if (nwords . le . 0) return 

Assign table indices 



PQQ 
PQG 
PGQ 
PGG 
PPL 
PMI 
PVA 



idPij(l,l) = 101 

idPij(2,l) = 201 

idPij(3,l) = 102 

idPij(4,l) = 202 

idPij(5,l) = 101 

idPij(6,l) = 101 

idPij(7,l) = 101 

Fill tables 

call MakeWRS(w, idPij(l,l) 

call MakeWtD(w, idPij(l,l) 

call MakeWtA(w, idPij(2,l) 

call MakeWtA(w, idPij(3,l) 

call MakeWtA(w, idPij(4,l) 

call MakeWRS(w, idPij(4,l) 

call MakeWtD(w, idPij(4,l) 

Done ! 

call clrUmsg ! clear s/r name for error messages 

return 

end 



PQQR, PQQS, AChi, 0) 

PQQD, AChi) 

PQGA, AChi) 

PGQA, AChi) 

PGGA, AChi) 

PGGR, PGGS, AChi, 0) 

PGGD, AChi) 



Figure 7: Subroutine that generates weight tables of user-given evolution kernels (LO only, in this 
example). The subroutine is passed to qcdnum via the routine f illwc, as is described in the text. 
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In Figure O we show the hsting of a custom weight routine (called myweight) that creates 
the weight tables of, in fact, un-polarised splitting functions in LO, see Appendix |X1 
The arguments of such a subroutine are as follows. 

subroutine mysub( w, nw, nwords, idpi j , mxord, idum ) 
implicit double precision (a-h,o-z) 
dimension w(*) , idPij(7,3) 



w The QCDNUM store (passed by reference); 

nw The number of words available in the store (input, passed by QCDNUm); 

nwords Number of words used by the tables (output); 

idpij List of weight table identifiers (output); 

mxord Maximum perturbative order supported by the tables (output); 

idum Not used at present. 

The body of the code in Figure [7] shows the steps to be taken in a custom weight routine. 

1. Set the maximum order of the custom evolution, here mxord = 1; 

2. Partition the store into tables by a call to booktab. Here are booked two type-1 
and two type-2 tables; 

3. Branch-out if there is not enough space in the store, as is signalled by a negative 
value of nwords returned by booktab; 

4. Return in idPij (id, lord) the identifiers of the various splitting function tables. 
The first index runs as follows 

12 3 4 5 6 7 



-T qq -< qg -< gq -< gg -< + -T- -TV 



There are only LO splitting functions in the example, with the non-singlet splitting 
functions all being equal to Pqq (table identifier 101); 

5. Generate the weight tables by calls to makewta, makewtb, makewrs, and makewd, 
as is described in Section 16.21 The calls in Figure [7] actually accommodate the 
splitting functions given in flA.3l) . 

One can have only one set of custom weight tables in memory; a second call to f illwc 
will result in a fatal error message. 

A custom evolution must obey the DGLAP evolution equations 02. 9p and 02.121) which 
implies that the evolution must properly split into singlet/gluon and non-singlet parts. 
We remind the reader that the evolution kernels in qcdnum are defined by convolution 
with parton number densities, and not parton momentum densities. These kernels can 
depend on x, Uf and /i^ and can be stored in tables of types-1, 2, 3 or 4. Note, however, 
that the /i^ dependence via as is taken care of in the evolution routine evolsg, through 
multiplication of the N^LO tables by (as/27r)^+^. 
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6.7 Error Messages in Add-On Packages 

When one writes an add-on package, the problem arises that QCDNUM error messages 
will be labelled with the name of the qcdnum routine and not with that of the package 
routine. This can of course become confusing for the user who should be aware only of 
the package routines, and not what is inside. 

One solution is that the package catches errors before qcdnum does, but this would 
duplicate a good checking mechanism which is already in place. An easier solution is to 
pass a string to qcdnum which contains the name of the package routine so that it will 
be printed together with the error message. For this, the routines setUmsg and clrUmsg 
are provided. For instance one of the first calls in the zmf illw routine of the zmstf 
package is 

call setUmsgCZMFILLWO 
so that, upon error, the user gets additional information: 



Error in BDDKTAB ( W, NW, ITYPES, NT, NWDS ) — > STOP 



No x-grid available 
Please call GXMAKE 

BOOKTAB was called by ZMFILLW 

The last call in zmf illw is 

call clrUmsg 

that wipes the additional message. This is important because downstream qcdnum 
errors would otherwise appear to have always come from zmf illw. 
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A Singularities 



In this appendix we denote by f{x) a parton momentum density and not a number 
density. In terms of / the convolution integrals in the evolution equations read 



/(x) 



[\zP{z)f(- 

J X ^^ 



>(0) 



The LO splitting matrix P^ in fl2.14p is written as, in the notation of (12. 7|) . 



p(0) p(0) 

-' qq -' qg 

p(0) p(0) 

-^ gq -* gg 



p(0) 9„ p(0) 

-Tqq ZUfT^q^g 

p(0) p(0) 

-* gqi -' gg 



(A.l) 



(A.2) 



The LO un-polarised splitting functions are given by 



P^l\^) 



^S(-) 



gg \-' 



4 
3 
1 

2 
4 

3 

6 



[x^ + (1 - x)2] 



l + (l-x 



X 



i2n 



X 1 — X I W Ti f 



(A.3) 



For the time-like evolution of fragmentation functions, the splitting functions Pq,g 
and Pgq, are exchanged in (]A.2p [15j. The '+' prescription in (1A.3P is defined by 



[f{x)U = fix)-6il-x) / fiz)dz 

Jo 



(A.4) 



so that 



For reference we give the expressions for /qq and Jgg obtained from (1A.3P and (lA.4p 



/(z)[^(^)]+d^= / [/(^)-/(l)](7(^)dz-/(l) / (?(^)d^. (A.5) 



i!Si^) 



4?( 



X 



3 ./,. 1 — 2; 



6 / dz 



X 



(l + z^)/(-j-2/(x)J+-/(x) 



+ 21n(l-x) 



6/(x) 



1-z 



ln(l - x) + 



^/ - -/(^) 



1 
+ Q I dz 

X 



1-z 



+ z(l-z) 



'X 



/&'- 



12 18 



(A.6) 



To write down a generic expression we decompose a splitting (or coefficient) function 
into a regular part (A), singular part (P), product of the two {RS) and a delta function 



P(x) = A(x) + [P(x)]+ + P(x)[^(x)]+ + K{x)5{l - x) 



(A.7) 
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where, of course, not all terms have to be present. The following functions are defined 
in the logarithmic scaling variable y = — ln(a;): 

hiy) = /(e^^), Qiy) = e-yp{e-y), A{y) = e'^ A{e-y) (A.8) 

with similar definitions for B and S] however, R[y) = -R(e~^) and K(y) = K{e^y) 
without a factor e~^ in front. With these definitions (lA.ip can be written as 

ry 
I{y) = / du Q{u) h{y - u) = h{y) + hiy) + hiy) + h{y) with 

ry 
hiy) = / du Aiu) hiy - u); 
Jo 

hiy) = r duBiu)[hiy-u)-hiy)]-hiy) r dzBiz); 
Jo Jo 

hiy) = f duSiu) [Riu)hiy-u)-RiQ)hiy)]-Ri^)hiy) [ dz Siz)- 
Jo Jo 

hiy) = Kiy)hiy). (A.9) 

where the last integrals of h and h are still expressed in the variable x = exp(— |/) to 
avoid integration extending to infinity in our expressions. Note that we are free to swap 
the arguments u and y — u in flA.Qp . 



B Triangular Systems in the DGLAP Evolution 

For the non-singlet evolution we have to solve the equation (see Section 13. 3p 

Va = b. (B.l) 

The matrix V is a lower triangular Toeplitz matrix, that is, a matrix with the ele- 
ments Vij depending only on the difference i — j as is shown in the 4x4 example fIS.lSp . 
This matrix is uniquely determined by storing the first column in a one-dimensional 
vector V so that Vij = fj-j+i for i > j, and zero otherwise. Eq. fIB.ip is, like any other 



lower triangular system, iteratively solved by forward substitution 

ai = h/vi 

i-l 

^^ for z > 2. (B.2) 



1 

Vl 



bi-^ V(^-j+l) a,j 



There is no recursion relation between aj_i and a^ so that in each iteration the sums must 
be accumulated, giving an operation count of nin + l)/2 for a system of n equations. 
This is as expensive (or cheap) as multiplying the triangular matrix by a vector. 

The substitution algorithm can be extended to solve the coupled singlet-gluon equation 



^m ^g& ) \9 ) \ ^ d, ) \ g ) V*'' 
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where a is a short-hand notation for Vqq, etc. These matrices are all lower triangular 
n X n Toeplitz matrices. Writing out this equation in components it is easy to see that 
for the first elements /i and gi we have to solve the 2x2 matrix equation 



(^i bi \ I fi \ _ I n \ ^ ( fi \ _ 1 I di -bi \ ( ri 



ci di J \ Qi J \ si J \ 91 J aidi - hci \ -Ci ai J \ Si 



(B.4) 



For i > 2 we have to accumulate the sums 



Ri = ^i-^^ ["(*+l-j) -fj + b(i+^'J) 9j] 

Si = Si-Y^ [c{i+i-j) fj + d(i+i-j) gj] (B.5) 



and solve, for each i, the equations 



ci di I \ Qi I \ Si 



(B.6) 



The operation count of this algorithm is four times that of (IB.2p . plus some little over- 
head to solve the 2x2 matrix equations for each i. 



C Zero Mass Structure Functions 

C.l General Formalism 

The zero-mass structure functions F2(x,Q^), Fi^{x,Q'^) and xF^.i^x.Q'^) in un-polarised 
deep inelastic scattering are calculated from (16.1 p with x = x. The Wilson coefficients 
are functions of a; (and sometimes Uf) only. We set, for the moment, the physical scale Q^ 
equal to the factorisation and renormalisation scale /i^ and write the singlet/gluon con- 
tribution to F2 and -Fl as (there is no contribution to xF^ since this structure function 
is a pure non-singlet) 

-j;(^^(x,g2) = [Q,s®gs](x,/i') + [C,,g®(7](x,/i2) ^ = 2,L. (C.l) 

Likewise, non-singlet contributions to the structure functions are given by 

-j;("^^(a;,Q2) = [a,ns®gns](x,/i') ^ = 2,L,3 (C.2) 

Jb 

where the label 'ns' stands for the non-singlet indices '+', '— ' and 'v' as defined by (12.111) . 
To be precise on notation: T2 = F2, T\^ = Fi^ and J-3 = xF^ in fIC.ip and flC.2l) . 



A structure function is calculated by adding the singlet/gluon and non-singlet parts, 
weighted by the appropriate combination of electroweak couplings; we refer to [28j for 
how to compute neutral and charged current cross sections and structure functions in 
deep inelastic charged lepton and neutrino scattering. 
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Like the splitting functions, the coefficient functions are expanded in powers of as, 

'^^'° = E«sCS^ ^ = 2,L,3 j = g,s,+,-,v (C.3) 



fc=0 



where I = (0, 1, 2) denotes (LO, NLO, NNLO) and a^ = aj2Ti. The LO coefficient 
functions are either zero or trivial delta functions: 



(C.4) 



cfl = 


<i = 5(1 - 


-x) 


cil = Sil - 


-x) 


^S = o 


cn = o 




ci% = 




^S = o 


ci:i = 




c!^L = 5(1 - 


-x). 



The NLO coefficient functions can be found in [H]. For those at NNLO we refer to 
[30| [311 [32] and the parametrisations given in [33] and [35] . 

The LO coefficient functions for Fl are zero so that the longitudinal structure function 
vanishes at LO. An alternative, which we call F{^, is calculated from the expansion 



fc=i 

In this way, Cl j is used already at LO (giving a non-zero Fl) and C\^ , at NLO. At 

NNLO the 3- loop coefficient function C^ j is taken from [55]. As stated in [35j, this 
3-loop calculation applies only to electromagnetic current exchange so that Z° or W^ 
contributions to F{^ at NNLO are, at present, not available. 



C.2 Renormalisation and Factorisation Scale Dependence 

To calculate the renormalisation scale dependence (/i^ ^ /i|) we replace, in the expan- 
sion of the coefficient functions, the powers of a^. by the Taylor series given in fl2.17p . If 
the expansion (JC.3P is used, the truncation of the right-hand side of (12.171) is to order a^ 



in NLO and a^ in NNLO. If, for F{^, the expansion (IC.SP is used, the truncation is to 
order a^ in LO, a^ in NLO and a^ in NNLO, like for the splitting functions. 

To calculate the factorisation scale dependence (Q^ 7^ Ai|), the coefficient functions 
in (JC.3P and (JC.51) are replaced by [33], [31] 



CfH ^ Cf^ and Cg^ ^ C^ + J^ Cg'"^^L^ fc > 1, (C.6) 



^(0) . r^(o) ^^A ni^) . r^W , '^ f^ik^Tn) 

m=l 

where Lp = \n{Q'^/fip) and /ip = /ip^. To write compact expressions for the ClJ , 
we introduce the following vector notation. In the non-singlet sector we have a one- 
dimensional vector Ci = Cj^ns and a 1 x 1 matrix P = P^s- In the singlet/gluon sector 
we have a 2-dimensional row- vector and a 2 x 2 matrix that are given by 



C, = (a,s a,g) and P- i^^^ ^ ^' 
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-* gq -* gi 



{k,m) 



In this vector notation, the functions C\-^ in (IC.6P are written as 






c 



(0) 

i 



>{0) 



C °) ® P(i) + Cf ^ ® P(°) - /3o / 



C 



(1,1) 



P(°) - /3o / 



C 



(3,3) 



cf ) ® p(2) + cf ) ® 


P 


« - A / 


+ cf^ ® 


p(0) 


-2/3o/ 


i{c!-^^ 


P(i) - /3i / 


+ Cf '^) ® 


P(°) - 2^0 J" 


} 


^cf ^^) « 


P(°) - 2/3o 


/ 













(C.7) 



t(2,2) 



For F2, Fl and xF^, the coefficients are calculated up to Cl^'^' . For F^, on the other 
hand, all coefficients in flC.7p are computed. Note, however, that quite some convolutions 
are trivial because the LO coefficient functions are either zero or (5-functions, see ( lC.4p . 

As mentioned above, the expression flC.6p applies only when /i| = /i^. It is therefore 
not possible to vary both scales /ip^ and Q^ at the same time. 



C.3 The ZMSTF Package 

The ZMSTF package is a qcdnum add-on with routines that calculate the structure 
functions F2, Fi^ and xF^ in un-polarised deep inelastic scattering. The structure func- 
tions are computed as a convolution of the parton densities with zero-mass coefficient 
functions, using the convolution engine described in Section [61 

The list of subroutines is given in Table El Note that error messages are, in most 
Table 6: Subroutine and function calls in ZMSTF. 



Subroutine or function 



Description 



ZMFILLW 


[ *nwords ) Fill weight tables 


ZMDUMPW 


[ lun, 'filename' ) Dump weight tables 


ZMREADW 


[ lun, 'filename', *nwords, *ierr ) Read weight tables 


ZMDEFQ2 


: a, b ) Define Q^ 


ZMABVAL 


[^ *a, *b ) Retrieve a and b coefficients 


ZMQFRMU 


[ qmu2 ) Convert /x| to Q^ 


ZMUFRMQ 


: Q2 ) Convert Q^ to /i| 


ZSWITCH 


[ iset ) Switch pdf set 


ZMSTFUN 


[ istf, def, X, Q2, *f , n, ichk ) Structure functions 



Output arguments are pre-fixed with an asterisk (*). 



cases, issued by the underlying qcdnum routines and not by the zmstf routine itself. 
However, the calling zmstf routine is mentioned in the error message so that you know 
where it came from. 
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call ZMFILLW ( *nwords ) 



Fill the weight tables. The tables are calculated for all flavours 3 < n/ < 6 and for all 
orders LO, NLO, NNLO. On exit, the number of words occupied by the store is returned 
in nwords. If you get an error message that the internal store is too small to contain the 
weight tables, you should increase the value of the parameter nzmstor in the include 
file zmstf . inc and recompile ZMSTF. 

This routine (or zmreadw below) should be called after an x-/x^ grid is defined in qcdnum 
and before the first call to zmstf un. 



call ZMDUMPW ( lun, 'filename' ) 



Dump the weights in memory via logical unit number lun to a disk file. The dump is 
unformatted so that the weight file cannot be exchanged across machines. 



call ZMREADW ( lun, 'filename', *nwords, *ierr ) 



Read weights from a disk file via logical unit number lun. On exit, nwords contains the 
number of words read into the store (fatal error if not enough space, see above) and the 
fiag ierr is set as follows. 

Weights are successfully read in. 

1 Read error or input file does not exist. 

2 Incompatible qcdnum version. 

3 Incompatible zmstf version. 

4 Incompatible x-/i^ grid definition. 

These errors will not generate a program abort so that one should check the value 
of ierr, and take the appropriate action if it is non-zero. 



call ZMDEFQ2 ( a, b ) 



Define the relation between the factorisation scale /ip and Q^ 

The Q^ scale can only be varied when the renormalisation and factorisation scales are 
set equal in qcdnum. The default setting is a = 1 and b = 0. The ranges are limited 
to . 1 < a < 10 and -100 < b < 100. 

A call to zmabval(a,b) reads the coefficients back from memory. To convert between 
the scales use: 

Q2 = zmqfrmu(qmu2) 
qmu2 = zmufrmq(Q2) 
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call ZSWITCH ( iset ) 



By default, the structure functions are calculated from the un-polarised parton densities, 
evolved with qcdnum (iset = 1). With this routine one can switch to the custom 
evolution (4), or to one of the external pdf sets (5-9). Switching to polarised pdfs (2) or 
to fragmentation functions (3) does not make sense and will produce an error message. 



call ZMSTFUN ( istf, def, x, Q2, *f , n, ichk ) 



Calculate a structure function for a linear combination of parton densities. 

istf Structure function index (1,2,3,4) = {Fi^,F2,xF3,F[). 

def (-6:6) Coefficients of the quark linear combination for which the structure func- 
tion is to be calculated. The indexing of def is given in (15. ip . 

x, Q2 Input arrays containing a list of x and Q^ (not fi"^) values. 

f Output array containing the list of structure functions. 

n Number of items in x, Q2 and f . 

ichk If set to zero, zmstf un will return a null value when x or /i^ are outside 

the grid boundaries; otherwise you will get a fatal error message. A /i^ 
point that is close or below the QCD scale A^ is considered to be outside 
the grid boundary. 

To calculate a structure function for more than one interpolation point, it is recom- 
mended to not execute zmstf un in a loop but to pass the entire list of interpolation 
points in a single call. The loop is then internally optimised for greater speed. 



D Heavy Quark Structure Functions 

A NLO calculation of the heavy quark contributions to the F2 and Fl structure functions 
in deep inelastic charged lepton-proton scattering is given in [16] . Only electromagnetic 
exchange contributions are taken into account. In this calculation, a heavy flavour h 
is not taken to be a constituent of the incoming proton but is, instead, assumed to be 
exclusively produced in the hard scattering process. Quarks with pole mass m < ruh are 
taken to be mass-less so that the input light quark densities should have been evolved 
in the FFNS with n/ = (3, 4, 5) for h = (c, b, t) [36] . 

A heavy flavour contribution to F2 or Fl is calculated from 

i^^^,Q^) = g{e^.^®C(>g(e^,®c2 + e^g.®C«+gp®Pi;))}, (D.l) 

where eh is the charge of the heavy quark (in units of the positron charge), g is the 
gluon density, qs is the singlet density and 

i=l 
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is the charge-weighted proton quark distribution for Uf hght flavours. The first term 
in flD.l|) is the LO contribution from the photon-gluon fusion process 'y*g — ?■ hh. The 



last three terms correspond to the NLO sub-process 'y*g — )■ hhg and 7*a — > hhqcj For 
the heavy quark coefficient functions C and V in fID.ip we refer to [16jo 

In terms of a number density /(x,/i^), the convolution integrals in (ID.1|) are defined by 

f®C= - zf{z, ^^') C{x/z, Q\ /i^ ml) (D.2) 

J ax ^ 

where a = l + 4m1/Q'^ and yU^ is the factorisation (equals renormalisation) scale which is 
usually set to /i^ = Q^ or /i^ = Q'^+Am\. The kinematic domain where the heavy quarks 
contribute is restricted by the requirement that the square of the 7*p centre of mass 
energy must be sufficient to produce the hh pair: W"^ = M^ + Q'^{l — x)/x > M^ -|- 4m^ 
so that the lower integration limit ax < 1 in (]D.2p . It turns out that the dependence of 
the coefficient functions on the relation between Q^ and /i^ cannot be factorised so that 
each setting of the scale parameters needs its own set of weight tables. To calculate the 
renormalisation scale dependence, the powers of Og = ol^JIt^ in (ID.ip are replaced by 
the Fourier expansion fl2.17p . truncated to Og in LO, and to Og in NLO. Note that one 
can vary either /i^ or Q^ with respect to //|, but not both at the same time. 

The convolution integral flD.2p is not of the general form (16.11) : (i) the factor x in front 
is missing; (ii) the pdf is xf{x) and not f{x) and (iii) the argument of C is x/z and 
not x/z- This mismatch is cured by presenting to QCDNUM the modified kernel 

C{X, A^^ Q^ TTT'h) = - C (-, Ai^ Q^ mi) , with X = ax. 
X Va / 

To make the heavy quark calculation available in qcdnum17 (as it was in qcdnumIG) 
we provide the add-on package hqstf described below. 



D.l The HQSTF Package 

The HQSTF package calculates up to NLO the heavy flavour contributions to the F2 or 
Fl structure functions from pdfs evolved in the ffns scheme with nj light flavours. The 
list of subroutines is given in Table [71 We will only describe here the routines hqf illw 
and hqstf un, the other ones being similar to those in the ZMSTF package. 



call HQFILLW ( istf, qmass, aq, bq, *nwords ) 



Fill the weight tables. To be called before anything else. 

istf Select structure function: 1 = Ft,, 2 = F> and 3 = both. 



^^In the LO and the first two NLO terms the virtual photon couples to the heavy quark, hence the 
factor el in (jD.l|) . The last NLO term describes the process where the virtual photon couples to a light 
quark which subsequently branches into a hh pair via an intermediate gluon; hence the appearance of 
the charge weighted sum, (jp, of light quark distributions. 

^^ Some of these coefhcient functions are given as interpolation tables (taken from code provided 
by S. Riemcrsma) since they are too complex to be cast into analytical form. Note that in 16 the 
coefficient functions are convolved with parton momentum densities and not with number densities. 
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Table 7: Subroutine and function calls in hqstf. 

Subroutine or function Description 

HQFILLW ( istf, qmass, aq, bq, *nwords ) Fill weight tables 

HQDUMPW ( lun, 'filename' ) Dump weight tables 

HQREADW ( lun, 'filename', *nw, *ierr ) Read weight tables 

HQPARMS ( *qmass, *aq, *bq ) Retrieve parameters 

HQQFRMU ( qmu2 ) Convert /i| to Q^ 

HQMUFRQ ( Q2 ) Convert Q^ to ii\ 

HSWITCH ( iset ) Switch pdf set 

HQSTFUN ( istf, icbt, def , x, Q2, *f, n, ichk ) Structure functions 
Output arguments are pre-fixcd with an asterisk (*). 



qmass (3) Input array with the c, b, t quark masses in GeV. If a quark mass is set to 

m/j < 1 GeV, no tables will be generated for that quark, 
aq, bq Defines the relation Q^ = afi^ + b. 

nwords Gives, on exit, the number of words used in the store. 

One will get a fatal error if the store is not large enough to hold all tables. In that 
case you can increase the value of nhqstor in the include file hqstf . inc and recompile 
HQSTF. The values of the mass and scale parameters can be retrieved at any time after 
the call to hqfillw (or hqreadw) by a call to hqparms(qmass,aq,bq). 



call HQSTFUN ( istf, icbt, def, x, Q2, *f, n, ichk ) 



Calculate the heavy quark contribution to a structure function. 

istf Calculate Fl (1) or F2 (2). 

icbt Select contribution from charm (l), bottom (2) or top (3). 

def (-6 : 6) Coefficients of the quark linear combination for which the structure func- 
tion is to be calculated. The indexing of def is given in f l6.8p . 

X, Q2 Input arrays containing a list of x and Q^ (not /i^) values. 

f Output array containing the list of structure functions. 

n Number of items in x, Q2 and f . 

ichk If set to zero, hqstf un will return a null value when x or /x^ are outside 

the grid boundaries; otherwise one will get a fatal error message. 

The routine checks that for icbt = (1,2,3) = (c,b,t) the pdfs were evolved in the ffns 
with Uf = (3, 4, 5) flavours. When icbt is pre-pended by a minus sign, the check on the 
FFNS remains active but that on the number of flavours is switched off. 

Here is a snippet of code that, in combination with zmstf, calculates the d,u,s contri- 
bution, the charm contribution and the total F2 (neglecting bottom and top) in charged 
lepton-proton scattering (the pdfs should have been evolved with n/ = 3 flavours). 
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dimension x(lOO) ,Q2(100) ,F2dus(100) ,F2c(100) ,F2p(100) 

dimension proton(-6:6) 

data proton /4.,1.,4.,1.,4,.1.,0.,1.,4.,1.,4.,1.,4./ ! divide by 9 

call zmstfun(2, proton, x, Q2, F2dus, 100, ichk) 
call hqstfun(2, 1, proton, x, Q2, F2c , 100, ichk) 
do i = 1,100 

F2p(i) = F2dus(i) + F2c(i) 
enddo 
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Index 



B-splines, 16-17 
backward evolution, 22-23, 31 
beta-functions, 6 
Bjorken-x variable, 8 

convolution integrals in QCDNUM, 18, 42 
convolution weights, 18 
custom evolution, 58-60 

DGLAP evolution equations, 7 
discontinuities in as evolution, 7, 13, 27 
discontinuities in pdf evolution, 13-14 

e^ basis pdfs, 11 

evolution of as, 6-7 

evolution of heavy quark pdfs, 14 

factorisation scale /ip, 8 

factorisation scale dependence, 65 

Fl structure function, 65 

FFNS fixed fiavour number scheme, 12 

flavour thresholds 

on the factorisation scale, 9 

on the renormalisation scale, 7, 9, 13 

forward substitution, 63-64 

Fourier convolution, 18 

fragmentation functions, see time-like evo- 
lution 

Gauss quadrature in QCDNUM, 44 

accuracy parameter of, 31 
generalised mass (gm) schemes, 42 

HQSTF package, 69-71 

indexing of e^ basis, 41 
indexing of q, q basis, 38 
interpolation mesh, 22 

ket notation, 10 

MBUTIL package, 24 

Mellin convolution, 8 

mesh point, see interpolation mesh 

MS scheme, 5 

multiple convolution, 18, 19, 47, 56, 57 

multiple evolution grid, 22, 27-29, 32 



non-singlet evolution, 9 
non-singlet quark density, 9, 10 
null value, 31 

number of active flavours nj, 6 
value at threshold, 9, 27, 52, 57 

parton density function, pdf, 6 

number/momentum density, 8 
parton luminosity, 19 
pdf type, pdf set, 35, 39, 40, 58 
PDG convention, 26 
PEGASUS program, 14, 27 
polarised evolution, 9 

QCDNUM program 

example job, 24-27 

execution speed, 29 

history, 5 

list of subroutines, 30, 44, 50, 54 

parameters in qcdnum.inc, 24 

program hang- up, 29 

web site, 24 
QCDNUM convolution engine 

booktab, 44 

copywgt, 47 

ef romqq, 51 

f crossf, 51 

f crossk, 50 

getalfn, 52 

getwpar, 49 

idspfun, 47 

makewrs, 46 

makewta, 45 

makewtb, 46 

makewtd, 46 

makewtx, 46 

nf lavor, 51 

qqfrome, 51 

scalewt, 47 

setwpar, 49 

stfunxq, 52 

tabdump, 49 

tabread, 49 

wcrossw, 47 

wtimesf , 48 
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QCDNUM fast convolutions 



f astclr, 54 


f astcpy, 58 


f astepm, 54 


f astfxf , 57 


f astfxk, 56 


f astfxq, 58 


f astini, 54 


f astkin, 57 


fastens, 55 


fastsum, 55 


QCDNUM routines 


asfunc, 37 


chkpdf 


40 


dmpwgt 


35 


evolfg 


38 


ffromr 


37 


fillwc 


58 


fillwt 


35 


fpdfij 


41 


fpdfxq 


41 


f snsij 


41 


f snsxq 


41 


f spine 


41 


f sumij 


41 


fsumxq 


41 


fvalij 


41 


fvalxq 


40 


getabr 


37 


getalf 


37 


getint 


31 


getord 


36 


getthr 


37 


getval 


31 


gqcopy 


34 


gqmake 


33 


grpars 


34 


gxcopy 


34 


gxmake 


32 


iqfrmq 


34 


ixfrmx 


33 


nwused 


36 


pdf inp 


39 


qcinit 


31 


qfrmiq 


34 


qqatiq 


34 


readwt 


36 



rfromf , 


37 


setabr, 


37 


setalf , 


37 


setint, 


31 


setlun, 


31 


setord, 


36 


setthr, 


37 


setval, 


31 


splchk, 


41 


xfrmix, 


33 


xxatix, 


33 



renormalisation scale fi'^, 6 
renormalisation scale dependence 

of pdfs, 10 

of structure functions, 65, 69 
rescaling variable X; 43 

scale parameter A, 7 

singlet quark density, 8 

singlet-gluon evolution, 8 

singlet /non-singlet basis e^, 11 

spline interpolation, 15-17, 41 

spline oscillation, 16, 22, 41 

splitting functions 

at leading order, 62 
perturbative expansion of, 9 
singularities in, 43, 62-63 
symmetries in, 8 

structure functions 

heavy flavour contributions, 68-69 
zero-mass structure functions, 64-66 

sum rule integrals, 26 

time-like evolution, 9, 62 

Toeplitz matrix, 19 

truncation prescription, 10, 52, 65 

t-variable, 18 

types of weight table, 43 

un-polarised evolution, 8 

VFNS variable flavour number scheme, 12 

weights, see convolution weights 
Wilson coefficient, 64 

y- variable, 18 

ZMSTF package, 66-68 
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